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ABSTRACT 


The  increase  in  fuel  prices  has  initiated  considerable 
interest  by  ship  operators  in  new  ship  autopilots  which 
minimize  the  propulsion  losses  due  to  steering. 

This  research  presents  the  results  of  work  on  a  steering 
design  for  the  Mariner  class  ship  based  on  a  computer  simu¬ 
lation.  A  model  of  the  Mariner  class  ship  was  coupled  to  a 
function  minimization  subroutine  to  minimize  the  added 
resistance  caused  by  rudder  activity  and  hull  drag  of  iner¬ 
tial  origins  caused  by  periodic  yawing  of  ship  in  seaway. 

The  Mariner  class  ship  computer  model  was  tested  in  calm 
water  and  in  a  seaway.  The  optimal  controller  parameters 
are  shown  in  look  up  tables  as  functions  of  ship  speed,  sea 
state,  encounter  angle  and  encounter  frequency.  This  tech¬ 
nique  can  be  used  as  an  adaptive  controller. 
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I.  INTRODUCTION 


Claims  by  many  researchers  indicate  that  a  carefully 
designed  controller  could  reduce  fuel  consumption  by  mini¬ 
mizing  the  propulsion  losses  which  are  caused  by  added  drag 
due  to  steering  of  the  ship. 

The  goal  of  this  thesis  was  aimed  at  developing  and 
demonstrating  the  utility  of  an  improved  steering  control 
for  the  Mariner  class  ship.  The  immediate  goal  of  this 
research  was  to  develop  design  methodology  for  an  adaptive 
autopilot  that  would  provide  effective  steering  control  with 
associated  cost  savings  for  a  full  range  of  seaway  and 
stability  conditions. 

To  simulate  the  ship  in  the  computer  program,  ship's 
nonlinear  equations  of  motion  were  needed.  Chapter  2 
addresses  the  Mariner  class  ship  nonlinear  equations  of 
motion. 

Chapter  3  addresses  the  work  on  testing  the  ship  simula¬ 
tion  model  and  open  loop  ship's  behaviors  in  calm  water  and 
in  a  seaway. 

The  basic  Nomoto  models  give  an  adequate  description  of 
ship  steering  dynamics  for  design.  The  Nomoto  second  and 
third  order  mo'1  2ls  were  developed  from  the  ship's  linear 
equations  of  mo  Chapter  4  adresses  the  Mariner  class 
ship  Nomoto  model  nations  by  using  mathematical  methods 
and  by  using  a  functio.  minimization  subroutine. 

Chapter  5  shows  an  adequate  cost  function  which  repre¬ 
sents  the  added  drag  due  to  steering  and  includes  deriva¬ 
tions  for  evaluating  the  weighting  factor.  Also  Chapter  5 
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X  =  A'-'X  +  B"U  then,  X(s) 
conditions  are  zero.) 


( s *1- A)'1  *B*U  (assuming  initial 


- 

V 

.  210 l”'s  + 

.0534 

R 

-  .  00  38 “S  - 

.0002 

D 


(7 . 67*s+ 1 )* ( 116 . 93*s+ 1 ) 


R(s)  =  s*YAW(s) 


so,  then 


YAW(s)  -0.189*(18.34*s+l) 


D ( s )  s*(7.67*s+l)*(116.93*s+l) 


result  is  K=0.189,  z=18.34,  Pl=7.67  and  P2=116.93 
Proceeding  to  the  second  order  Nomoto  equation: 
YAW ( s  )  K 


D ( s )  s*(Pl*s  >  1) 

Deriving  the  second  order  Nomoto  transfer  function  from 
the  yaw  equation  only,  the  result  is  K=0.03  and  Pl=10. 

B.  COMPUTER  APPROACH 

We  used  a  function  minimization  subroutine  to  obtain 
parameters  of  the  transfer  functions.  Figure  4.1  shows  the 
scheme  used  to  obtain  the  third  and  second  order  Nomoto 
transfer  functions.  The  results  are  given  in  Table  IV. 
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IV.  NOMOTO  MODEL  OF  THE  MARINER 


To  find  a  model  which  can  be  used  in  computer  simulation 
the  Mariner's  linear  equations  of  motion  and  its  hydrody¬ 
namic  coefficients  were  used  and  third  and  second  order 
Nomoto  transfer  functions  derived.  Values  used  were  for  ship 
speed  of  15  knots. 

Mariner’s  linear  equations  of  motion  are 

(m  -  Xa)*U  =  X u" DU  ( eqn  4.1) 

(m  -  Y*)*V  -  Y v*V  =  (Y—  m*X6)*R  ♦  (YR-  m*Ul)*R 
(I  -  Nj,  )*R  -  (Nr-  m*XG*Ul)*R  =  (N*-  m*XG)*V  +  Nv*V 

A.  MATHEMATICAL  APPROACH 

Proceeding  to  the  third  order  Nomoto  equation: 

YAW ( s )  K* (Z*s  ♦  1) 


D(s )  s*(Pl*s  ♦  l)*(P2*s  ♦  1) 

Deriving  the  third  order  Nomoto  transfer  function  from 
the  sway  and  yaw  equations,  we  show  them  in  matrix  notation 
as  follows. 


V 

_ 

- .0372 

-8.42 

... 

V 

+ 

.210 

R 

- .0003 

-  .  10 

R 

-  .003 
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The  disturbance  forces,  moments,  and  added  mass,  added 
inertia  terms  were  found  by  running  the  sea  state  program 
that  is  presented  in  [Ref.  4].  Data  for  a  Mariner  class 
ship  and  chosen  sea  conditions  that  were  used  in  the  sea 
state  program  are  shown  in  APPENDIX  B. 

For  observation  purposes  'FX'  (disturbance  force  in 
surge),  'FY'  (disturbance  force  in  sway)  and  'MZ'  (distur¬ 
bance  moment  in  yaw)  were  added  into  the  surge,  sway  and  yaw 
equations  that  were  used  in  the  simulation  program. 

For  regular  seas  the  program  has  been  run  a  few 
times. Every  time  different  FX,  FY,  MZ  and  coefficients  were 
used  to  represent  different  ship  characteristics  such  as 
ship  speed,  loading  2  etc.  and  enviromentel  conditions  such 
as  sea  state,  encounter  angle,  encounter  frequency.  Results 
show  that  U,  V  and  R  are  sine  waves  with  amplitude  and  phase 
depending  on  ship  and  environmental  conditions. 

Time  response  of  U,  V  and  R  in  200  seconds  are  presented 
in  Figure  (  3.2),  for  ship  speed  15  knots,  encounter  angle 
030.0  degree,  encounter  frequency  0.60  radian/ second  and  sea 
state  6. 


2During  this  research,  displacements  (Molded)  up  to  32 
ft  were  chosen  as  a  loading  condition  for  the  Mariner. 
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III.  OPEN  LOOP  SHIP’S  BEHAVIOR  IN  CALM  WATER  AND  SOME  SEA 

STATES 


A.  CALM  WATER  CASE 

Using  a  ship's  nonlinear  equations  of  motion  and  Mariner 
class  ship  coefficients,  a  simulation  program  THESIS  FORTRAN 
was  developed  and  run  to  observe  U,  V  and  R.The  computer 
program  THESIS  FORTRAN  is  shown  in  Appendix  A. 

First  run  D=0,  YAWC=0,  Ul=15  knots  were  applied  and  it 
was  seen  that  the  ship  stays  on  its  initial  course  and 
speed.  U=15  knots,  V=R=0. 

The  program  was  rerun  a  few  times,  changing  the  rudder 
angle  to  2.8  degrees  to  both  sides  (port  and  starboard)  and 
it  was  observed  that  increasing  the  rudder  angle  changes  the 
ship's  course  and  also  U  decreased,  the  absolute  values  of  V 
and  R  increased.  After  a  few  hundred  seconds  U,  V  and  R 
reached  steady  values  independently.  These  steady  values 
depend  on  rudder  angle,  large  rudder  angles  decrease  U  and 
increase  V  and  R.  For  a  rudder  angle  of  one  degree  and 
constant  speed  (Ul)  of  15  knots,  the  first  200  seconds  of 
time  response  of  U,  V  and  R  are  shown  in  figure  3.1  as  an 
example . 

B.  SEA  STATE  CASE 

To  observe  the  behavior  of  the  ship  in  a  sea  state, 
disturbance  forces  and  moments  are  needed  that  depend  on  sea 
state,  ship  speed,  encounter  angle  and  encounter  frequency. 
Also  in  sea  steate  hydrodynamic  parameters  are  changed,  i.e, 
the  added  mass  and  added  inertia  are  functions  of  encounter 
frequency  and  sea  state. 


20 


TABLE  III 

Assessment  of  the  Coefficients  in  the  N_Equation 


Variable 

Coefficient 

Value  of  Coeff: 

• 

V 

(m*X6-  N-) 

-0.00478 

• 

R 

(IX‘ 

0.0175 

V 

N. 

-0.0555 

V3 

(l/6)*Nyvv 

0.345 

V*R2 

0.5*N^k^ 

0 

V*D2 

0.00264 

V*DU 

Nvvx 

_ 

V '-DU  2 

0  •  5  -  N  v 

R 

(Nr-  m*X6*Ul) 

-0.0349 

R3 

(1/6)"N^^ 

0 

R"  V  2 

0 . 5 

-1.158 

R*  D  2 

°-5*N^X> 

0 

R-'DU  • 

N*vx 

_ 

R  "DU  2 

0 . 5 

_ 

D 

Nt, 

-0.0293 

D3 

(l/6)*»m 

0.00482 

D*V 2 

0.5*Nt>vv 

0.1032 

D*R2 

°-5*Norr 

0 

D--DU 

0 

D  "DU 2 

0-5**W 

_ 

V*R*D 

0 

N° 

0.00059 

DU 

Na 

0 

DU2 

N° 

11  SAVA. 
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TABLE  II 


Assessment  of  the  Coefficients  in  the  Y_Equation 
Variable  Coefficient  Value  of  Coefficient 


V 

(m  -  Y<) 

0.327 

R 

(m*XG-  Y-  ) 

-0.0018 

V 

-0.244 

V1 

d/6)^Vvv 

-1.702 

V“R2 

°*5*Yv^ 

0 

V*D  2 

0-5*Yv>t> 

-0.0008 

V'-DU 

Y 

L  NiU. 

V*DU 1 

°-5"Yy^ 

R 

(Yk-  m) 

-0.105 

R1 

0 

R*V  2 

°-5*Yrvv 

3.23 

R*D  2 

°-5*Ym 

0 

R*DU 

Y^ 

R"  DU  2 

°-5*Y*vxv* 

__ 

D 

Y0 

0.0586 

D 1 

(l/6)*Yeoo 

-0.00975 

D*V  2 

°-5*YBvv 

0.25 

D*R2 

0-5“Yo^ 

0 

D*DU 

Yt,vx 

0 

D*DU  2 

V*R*D 

°‘5*Yt>vxvx 

YyR.0 

Y® 

0 

-0.0008 

V*R*D 


0 

-0.0008 


TABLE  I 

Assessment  of  the  Coefficients1  in  the  X_Equation 


Variable 

Coefficient 

Value  of  Coeff: 

• 

U 

m  - 

0.177 

DU 

Xvx 

-0.0253 

DU2 

0.5*X^ 

0.00948 

DU3 

(1/6)“X^^ 

-0.00217 

V2 

0.5*XWV 

-0.189 

R2 

0.5*Xt?>+  m*X& 

0.00379 

D2 

0 . 5*X0X) 

-0.02 

V 2  *DU 

0 . 5  ••  X  v  v  ^ 

— 

R 2  *DU 

— 

D  2  *DU 

0.5"XOOVX 

— 

V*R 

X  m 

0.168 

V*D 

xv0 

0.0196 

R*D 

0 

V*R*DU 

X\/RW 

— 

V*D*DU 

XVDMk 

— 

R*D*DU 

X  RSv\ 

_ 

x° 

0 

lAll  derivatives  are  nondimensionalized  on  the  basis  of 
RHO , L , T  and  S. 

No  entry  in  these  columns  means  the  coefficient  was  ignored. 


These  surge,  sway  and  yaw  equations  can  be  rewritten  in  the 
form: 


dU/dt  =  gt(t,  U  ( t )  ,  V(t),  R  ( t )  ,  D(  t ) ) 


dV/dt  =  gj(t,  U  ( t ) ,  V  ( t )  ,  R(t),  D  ( t ) ) 


(eqn  2.10) 


dR/dt  =  g^(t,  U(t ) ,  V(t ) ,  R(t),  D  ( t ) ) 

Where  U(t),  V(t),  R(t)  and  D(t)  are  the  instantaneous  values 
of  U,  V,  R  and  D  at  any  time  t. 

Equation  2.10  is  a  set  of  three  first-order  differential 
equations  for  which  approximate  numerical  solutions  are 
readily  obtained  on  a  digital  computer.  The  key  to  the 
numerical  solution  is  that  values  of  U,  V  and  R  at  time 
t+DELT  are  obtained  from  knowledge  of  the  values  of  U,  V,  R 
and  D  at  time  t  using  a  simple  first-order  expansion;  that 
is, 

U(t+DELT )  =  U(t)  +  DELT*U(t) 


V (t+DELT )  =  V(t)  +  DELT *V ( t ) 


(eqn  2.11) 


R( t+DELT )  =  R(t )  +  DELT*R(t) 

This  method  is  found  to  give  adequate  accuracy  for  the 
present  type  of  differential  equations  because  of  the  fact 
that  the  accelerations  U,  V  and  R  vary  but  slowly  with  time, 
due  to  the  large  mass  and  inertia  of  a  ship  compared  to  the 
relatively  small  forces  and  moments  produced  by  its  control 
surface.  Any  desired  accuracy  of  the  solutions  can  be 
obtained  with  a  computer  by  simply  using  smaller  time  inter¬ 
vals  DELT.  This  procedure  was  used  for  all  computer  programs 
which  were  developed  for  this  thesis. 


V=D=0.0  are  identified  as  Y° and  N*  these  are  likely  to  be 
speed  dependent.  To  see  the  rudder  effects  in  calm  water  or 
sea  state  propeller  effects  were  ignored. 

Finally,  from  the  X  ,  Y  and  N  equations  the  ship's  surge, 
sway  and  yaw  equations  can  be  written  as  follows. 


fjfU.V.R.D) 


(eqn  2.7) 


(m  -  XJ 


(Iz-  Nk)*f2(U,V,R,D) 


(m*X6-  Yk)*f^U,V,R,D) 


(eqn  2.8) 


(m  -  Y.)*(IZ-  N^)  -  (m*Xc-  Y • ) 


Y.  )*f5(U,V,R,D) 


-  (m *X&-  N;)*f(U,V,R,D) 


(eqn  2.9) 


Where : 


Xuu=  d2X/dU2  X6„  =  d2X/dD2  X^=  d2X/dU*dV 

X^  d 3X/ dU 3  Xuu=  d3X/dU2*dR  etc. 


f  U ,  V ,  R ,  D )  =  0 . 5*XWJk*DU2  ♦  (l/6)*XxWkU*DU*  (eqn  2.4) 

♦  0 . 5  *X  vv  *V 2  +  (0.5*XKjL+  m*X6)*R2 
+  0.5*Xvvu.*V2*DU  +  0.5*X1l1Uj*R2*DU 
+  m)*V*R  ♦  Xvo"'V"'D  ♦  Xvll>*V*R*DU 

+  XVOvJ*V*D*DU  +  X(l^*R*D*DU  +  X^VDU  ♦  X° 

+  0.5*XW*D*  -  0 . 5 “X 2 *DU  +  X„*R*D 


f^U.V.R.D)  =  "DU  ♦  Y^*DU'+  Yv*V  +  Y0*D  (eqn  2.5) 
+  Y«  +  0.5*YvM*V*R2  -  0 . 5 *Y ^ 0 ^ * V * D 2 
+  Y v  *V -DU  ♦  0.5*Yv^*V*DU2  +  (Y--  m*Ul)*R 


>  (l/6)*Ym*D3  ♦  Y^*R*DU  ♦  0.5*Y^*R*DU* 

+  (l/6)*YRiR*R*  +  0.5*Y*v>(*D*Va  ♦  0 . 5*Yg.vv*R*Vi 
+  0-5*Yd|1r*D*R2  ♦  Y0**D*DU  -  0.5*Y^*D*DU2 
+  YwR0*V*R*D  +  (l/6)*Yvyv*VJ  *  0.5*Ymd*R*D‘ 


f  (U  ,  V  ,  R ,  D )  =  N0^  -'DU  +  N°USX*DU2  +  Nv*V  ♦  Nft*D  (eqn  2.6) 

♦  0.5*NPNk^*D*DU2  +  ( Ng - m*X<- U 1 ) *R 

♦  NVv*V*DU  +  0 . 5*N ^vjL\jfrV*DU 2  ♦  NvjL9*V*R*D  ♦  N° 

♦  0.5"N0VV"D"V2  +  0.5*N6Rh,*D*R2  *  Nfty*D*DU 

♦  NRvk*R*DU  ♦  0.5*NRw*R*DUa  ♦  (l/6)*NDfiS*D* 

+  (1/6)*N^*RJ  +  0.5*NRvV*R*V2  ♦  0.5*Nrw*R*D2 

♦  (l/6)*Nvvv*V3  +  0 . 5*NvRR*V*R2  ♦  0.5*NVOI*V*D: 


All  of  the  derivative  coefficients  of  the  equations  are 
evaluated  on  the  basis  of  experimental  data  obtained  from 
captive  model  tests,  and  given  in  Table  I,  II  and  III, 

The  Y  force  and  N  moment  induced  by  the  rotation  of  a 
single  propeller  or  by  unirotating  multiple  propellers  at 


II.  NONLINEAR  EQUATIONS  OF  MOTION 


Nonlinear  equations  of  motion  are  suitable  for  predicting 
tight  maneuvers  and  also  suitable  for  computer  programming. 
The  nonlinear  equations  of  motion  used  i.i  this  work  have 
been  developed  by  Abkowitz  [Ref.  1,  2],  and  Strom_Tejsen 
[Ref.  3],  based  on  a  Taylor  series  expansion  of  forces  and 
moments.  Terms  higher  than  third  order  are  not  included  in 
the  equations  because  experience  has  shown  that  accuracy  is 
not  significantly  improved  by  their  inclusion. 

A  result  of  symmetry  about  the  xz-plane,  X  is  an  even 

•  • 

function  of  V,  R,  D,  V  and  R  so  on,  the  crosscoupled  terms 
in  the  equations  involving  odd  powers  of  V,  R  and  D  are 
zero,  however,  crosscouple  terms  which  involve  even  powers 
of  V , R  and  D  are  nonzero.  In  contrast  to  X,  the  expressions 
for  Y  and  N  are  odd  functions  of  V,R,D,V  and  R;  that  is, 
only  the  coefficients  of  the  terms  in  the  expansion  with  odd 
powers  are  nonzero;  those  with  even  powers  are  zero.  For 
some  reasons,  X  is  neither  an  odd  nor  an  even  function  of  U 
but  rather  its  expansion  includes  all  powers  of  DU. 

Equations  X,Y  and  N  are  functions  of  U,V,R,U,V,R  and  D. 
Taylor  series  expansions  of  X,Y  and  N  including  terms  up  to 
the  third  order  are  as  follows: 


(m  -  X*)*U  =  fj( U , V , R , D ) 

(m  -  Y*)*V  ♦  (m*Xc-  Yk)*R  =  f/U,V,R,D) 
(m*Xfi-  N*)*V  ♦  (I2-  Nk)*R  =  f U  ,  V  ,  R , D ) 


(eqn  2.1) 
(eqn  2.2) 
(eqn  2.3) 
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presents  the  assumptions  and  approaches  needed  to  find  the 
cost  function  that  is  used  by  many  researchers. 

Ship  dynamics  change  with  operating  conditions  such  as 
ship  speed,  encounter  angle  and  encounter  frequency. 
Chapter  6  presents  optimal  controller  parameters  as  a  func¬ 
tion  of  different  operating  conditions. 

Chapter  7  addresses  an  approach  to  an  adaptive  controller 
utilizing  information  which  is  easy  to  measure  on  ship  board 
such  as  ship  speed,  heading  error  and  rudder  angle.  This 
adaptive  controller  must  be  used  to  provide  minimum  added 
drag  due  to  steering. 

Conclusions  were  drawn  from  simulation  results  and  are 
presented  in  Chapter  8.  This  chapter  also  recommends  some 
future  studies,  which  can  be  done  as  extensions  of  this 


TABLE  IV 

The  Nomoto  Model  Parameters  for  Mariner 
Third  Order  Second  Order 

K  =  0.189  K  =0.0298 

Z  =  18.347  PI  =  9.989 

PI  =  7.6739 

P2  =  116.929 

As  seen  the  answers  obtained  by  function  minimization 
agree  closely  with  the  analytic  solutions. 
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V.  DERIVATION  OF  A  COST  FUNCTION  FOR  THE  MARINER  SHIP 


In  recent  years  many  researchers  have  studied  the  problem 
of  optimizing  an  automatic  ship  steering  controller  for 
minimum  fuel  consumption.  It  is  well  known  that  additional 
drag  is  introduced  by  steering  and  that  both  the  rudder 
motion  and  the  yawing  motion  contribute  to  this  added  drag. 

When  deriving  a  cost  function  we  are  required  to  find  one 
cost  function  that  must  be  convenient  for  ship  board  use. The 
cost  function  that  is  commonly  used  in  recent  years  is 


J  =  /  (LAMDA'VYAWE2  +  D2  )*dt 

J 

o 

Where  YAWE  =  Yaw  error 

LAMBDA  =  Weighting  factor 


(eqn  5.1) 


Derivation  of  this  cost  function  from  the  surge  equation  has 
been  well  explained  in  [Ref.  5]  for  the  SL/7  class  ship  by 
R.E.REID. 


To  derive  the  cost  function  for  the  Mariner,  REID's 
approach  is  taken  as  a  reference  and  his  assumptions  are 
used.  To  show  how  the  cost  function  for  the  Mariner  was 
derived,  REID's  work  is  presented  here  step  by  step  with 
derivations  of  the  Mariner's  cost  function. 

The  Surge  Equations: 

SL/7  : 


(m  -  XJ*U  =  X°  +  0 . 5*Xvv*V*  ♦  0.5*Xoft*D*  (eqn  5.2) 

+  (Xva>  +  m)*V*R  ♦  Xp 

Where:  X°  =  -0.0003  for  SL/7 

Xp  =  Propeller  thrust 
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MARINER: 


(m  -  Xa)*U  =  0 . 5  *X  vv  *V 2  +  (XVR>  ♦  m)*V*R 

+  0 . 5  *X  *D  2  +  ( 0 . 5  *X  ^  ♦  m*X  6 ) *R s 
+  Xv0*V*D  +  X^'-'DU  ♦  0.5*Xwk*DU2 
+  ( 1  /  6  )  *X  vmjjoT  DU 1 


(eqn  5.3) 


It  is  seen  that  there  are  some  terms  in  the  Mariner's 
surge  equation  which  they  are  not  included  in  the  SL/7's 
surge  equation.  Assuming  steady  state  situations  since,  U=0. 
The  instantaneous  surge  relevant  to  steering  is 


DX  =  0.5*XVn*V2  +  0 . 5*X^0*D2  +  (Xv^  +  m)*V*R  (eqn  5.4) 

SL/7: 

MARINER: 

DX  =  0 . 5  *X  *V 2  *  0.5*XW*D*  ♦  (XVB>  +  m)*V*R  (eqn  5.5) 
+  (0.5*Xk?l  +  m*X  ) *R 2  ♦  Xm0*V*D  +  X^-DU 
+  0.5*X^*DU2  +  ( 1  /  6  )  *X  v^vxvxVDU 3 


From  the  instantaneous  surge  equation  relevant  to 
steering  of  the  SL/7  Reid  came  up  with  the  following  cost 
function: 


J  =  (1/2T)*  / (LAMBDA" *V*R  ♦  ETA*V 2  ♦  D2)*dt 


(eqn  5.6) 


o 

Where 


LAMBDA" = 

(m 

*  Xvr)/0.5*Xw 

ETA 

(0. 

5*Xvv)/0.5*Xm, 
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and  he  found  that  the  value  of  the  (ETA-’-'V2)  term  is  very 
small  so  he  neglected  the  (ETA-’-V2)  term,  then  the  cost  func¬ 
tion  for  the  SL/7  is 


■dt  (eqn  5.7) 

Using  the  same  approach  for  the  Mariner,  the  following 


cost  function  was  derived. 


J  =  (1/2T)*/(A1*DU  +  A2--DU2  +  A3  “DU 3  +  A4--V2  (eqn  5.8) 


-  A5 *R 2  +  D2  +  A7*V*R  -  A8*V*D)*dt 


Where  Al  = 

Xu/0.5*Xftt> 

A2  = 

0. 

,5*X^/0.5*X> 

A3  = 

(1/6)*X^/0.5*Xw 

A4  = 

0. 

.5*Xvv/0.5*Xy 

A5  = 

(0.5*X?*  *  m*XG)/0. 

A7  =  (Xv*  *  m)/0.5*Xt*  A8  =  X>,o/0.5*Xw 

A5  and  A8  are  always  negative  numbers,  since  then  A5  and  A8 
have  a  minus  sign  in  the  equation.  For  the  calm  water  case 
and  when  Ul=15  knots,  D=2.6  degrees  after  2000  seconds, 
values  of  every  term  in  the  Mariner's  cost  function  are 
given  below  to  give  an  idea  about  assumptions. 


A1*DU  = 

-0.001939 

A  2  "'DU  2  = 

-0.00000111 

A3*DU3= 

-0.00000000039 

A4“V2  = 

0.000003253 

A5*R2  = 

-0.0002884 

A7*V*R  = 

0.0001923 

D2 

0.002059 

n 

D 

> 

* 

OO 

C 

-0.00002609 

As  seen  from  the  above  (A2*DU2),  (A3*DUJ),  (A4*V2)  and 

( A8 “V"D )  terms  are  very  small  compared  to  others,  so  they 


may  be  neglected.  Also  to  measure  the  DU  on  shipboard  is 
very  hard  although  it  may  someday  be  measured  by  new  satel¬ 
lite  facilities,  so  we  must  not  include  it  in  the  cost  func¬ 
tion.  After  these  assumptions  the  cost  function  for  the 
Mariner  is 

rK 

J  =  0.5*I(-A5*R2  +  D2  +  A7*V*R)*dt  (eqn  5.9) 
o 

The  only  difference  between  Equation  (5.9)  and  Equation 
(5.7)  is  the  (A5'"R2 )  term  that  is  not  included  in  the  SL/7's 
cost  function.  To  see  the  effect  of  the  (A5*R2 )  term  on  the 
cost  function  the  Mariner's  cost  function  was  evaluated  with 
and  without  the  (A5"“R2  )  term  for  the  calm  water  case  and 
Ul=15  knots,  D=2.6  degree  after  2000  seconds,  results  are: 


J 


with  ( A5*R2 ) 

0.0025399826296 

without  (A5,VR2) 

0.00 22515066462 

There  is  no  big  difference  between  these  two  J  values,  and 
to  make  the  derivations  similar  to  Reid's  derivations  the 
(A5‘VR2 )  term  won't  be  included  in  the  Mariner's  cost  func¬ 
tion  but  as  it  is  known  that  for  the  Mariner  the  (A5*R2) 
term  is  as  big  as  the  j(A7*V*R)  term,  it  would  be  better  to 
consider  it  in  the  cost  function.  After  all  of  the  above 
steps  the  cost  function  of  the  Mariner  may  be  written  as  in 
Equation  (5.7). 


V  and  R  are  hard  to  measure  on  shipboard,  but  (V"R)  can 
be  defined  as 

V*R  =  OP "YAW1 

Where  R  =  YAW  =  YAWE*w 

OP  =  Distance  from  the  ship  pivot  point  to 
the  origin.  =  0.3*L 

w  =  Natural  frequency  of  the  ship's  steering 
system  closed  loop,  (w  =  0.05  rad/sec  was  initially  used.) 

Finally  the  cost  function  for  the  Mariner  is 


J  =  0.5* /(LAMBDA* YAWE 2  +  D2 )*dt 


( eqn  5.10) 


Where 


LAMBDA  =  (m  +  X^  )*0P*w*  /  0 . 5*XW 


Since,  X  depends  on  ship  speed,  for  different  ship  speeds 
values  of  LAMBDA  were  calculated  and  are  presented  in  Table 
(V).  These  LAMBDA  values  were  calculated  by  assuming  the 
natural  frequency  of  the  ship's  steering  system  closed  loop 
is  equal  to  0.05  radian/ second .  How  important  the  accuracy 
of  the  LAMBDA  value  is  with  respect  to  finding  the  optimal 
control  parameters  will  be  observed  in  Chapter  6. 


TABLE  V 


Values  of  Weighting  Factor  for  Different  Speeds 
of  The  Mariner  Class  Ship 


ship  speed  | 

LAMBDA 

(Knots)  | 

10  | 

6.57 

15  I 

2.91 

20  | 

1.64 

VI.  CONTROLLER  DESIGN  FOR  THE  MARINER,  USING  FORTRAN  PROGRAM 


We  coupled  a  function  minimization  subroutine  to  the 
cascade  compensater  that  is  coupled  with  Mariner's  equations 
of  motion  and  used  the  subroutine  to  adjust  the  controller 
parameters  to  minimize  the  cost  function  which  is  derived  in 
Chapter  V,  and  evaluate  the  minimum  cost. The  Fortran  program 
to  calculate  the  optimal  parameters  is  given  in  Appendix  C. 
Compensater  'A'  and  'B'  are  used  as  controllers,  their 
structures  are  shown  in  Figure  (  6.1  ). 


Kx(SxZ1  +  1) 

(SxPI  f  1) 


COMPF  N^AT  OR  -  A 


Figure  6.1  Various  Structure  Controllers. 


A.  CALM  WATER  CASE 

For  calm  water,  for  a  given  1  degree  yaw  command,  using 
the  computer  method  we  optimized  controllers  'A'  and  'B'  and 
the  results  are  shown  in  Tables  VI,  VII  and  VIII 


TABLE  VI 

Optimal  Controller  Parameters  I 


Simulation  Results  -  Steady  State  600  seconds 

For  ship  speed  10  Knots  the  optimal 
parameters  of  various  controllers 


and  the 

cost 

CONTR 

K 

Z1 

PI 

P2 

COST 

A 

0.77 

60.07 

20.46 

0.0494896 

B 

0.74 

60.07 

20.12 

0.902 

0.0514841 

TABLE  VII 

Optimal  Controller  Parameters  II 

Simulation  Results  -  Steady  State  600  seconds 

For  ship  speed  15  Knots  the  optimal 
parameters  of  various  controllers 


and  the 

cost 

CONTR 

K 

Z1 

PI 

P2 

COST 

A 

0.53 

51.46 

18.08 

0.0186974 

B 

0.50 

51.58 

17.81 

0.890 

0.01957901 

TABLE  VIII 

Optimal  Controller  Parameters  III 

Simulation  Results  -  Steady  State  600  seconds 

For  ship  speed  20  Knots  the  optimal 
parameters  of  various  controllers 


and  the  cost 

CONTR  K  Z1  PI  P2  COST 

A  0.40  44.87  16.06  -  0.0094217 


B  0.39  41.11  15.84  0.880  0.00991801 


These  results  will  be  references  for  the 
design  for  sea  state  operation.  We  observe  that 
the  speed  gives  us  smaller  controller  parameters, 
behavior  can  be  seen  from  the  SL/7's  results 
presented  in  [Ref.  6]. 


controller 
increas ing 
also  this 
that  are 


B.  REGULAR  SEAS  CASE 

The  ship  in  regular  seas  is  affected  by  sea  wave  distur¬ 
bance  forces  and  moments.  These  are  functions  of  sea  state 
and  encounter  frequency.  Also  the  added  mass  and  the  added 
inertia  terms  are  functions  of  sea  state  and  encounter 
frequency,  and  the  encounter  frequency  depends  on  the 
encounter  angle  and  ship  speed.  All  of  these  variables  must 
be  considered  when  calculating  the  optimal  parameters  of  the 
controller . 

Using  the  computer  method  controllers  'A'  and  'B'  were 
optimized  for  a  few  different  cases  and  the  results  are 
shown  in  Tables  IX,  X,  XI,  XII,  XIII,  XIV,  XV,  XVI  and  XVII. 


TABLE  IX 

Optimal  Controller  Parameters  IV 

Simulation  Results  -  Steady  State  600  seconds 
optimal  parameters  of  various  controllers  and  the  cost. 

For  ship  speed  15  Knots,  encounter  angle  030.0  degree, 
encounter  frequency  0.50  rad. /sec.  ana  sea  state  6. 


CONTR 

K 

Z1 

PI 

P2 

COST 

A 

0.358 

6  6.6 

24.61 

0.001252939 

B 

0.35 

44.68 

06 . 58 

9  .  720 

0.0010086581 
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TABLE  X 

Optimal  Controller  Parameters  V 


Simulation  Results  -  Steady  State  600  seconds 
optimal  parameters  of  various  controllers  and  the  cost. 

For  ship  speed  15  Knots,  encounter  angle  060.0  degree, 
encounter  frequency  2.50  rad. /sec.  ana  sea  state  9. 


CONTR  K  Z1 

A  0.60  41.58 

B  0.70  30.49 


PI 

P2 

COST 

11.33 

0.000020747 

03.37 

2.940 

0.000016038 

TABLE  XI 

Optimal  controller  Parameters  VI 


Simulation  Results  -  Steady  State  600  seconds 
optimal  parameters  of  various  controllers  and  the  cost. 

For  ship  speed  15  Knots,  encounter  angle  090.0  degree, 
encounter  frequency  0.50  rad. /sec.  ana  sea  state  7. 


CONTR  K  Z1  PI  P2  COST 

A  0.54  30.65  09.35  -  0.054223988 


B  **********  IT  DID  Not  CONVERGE  ************ 


TABLE  XII 

Optimal  Controller  Parameters  VII 


Simulation  Results  -  Steady  State  600  seconds 
optimal  parameters  of  various  controllers  and  the  cost. 

For  ship  speed  15  Knots,  encounter  angle  120.0  degree, 
encounter  frequency  0.75  rad. /sec.  and  sea  state  7. 


CONTR 

K 

Z1 

PI 

P2 

COST 

A 

0.52 

41.09 

08.95 

0.018345 

B  **********  IT  DID  NOT  CONVERGE  ************ 
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TABLE  XIII 

Optimal  Controller  Parameters  VIII 

Simulation  Results  -  Steady  State  600  seconds 
optimal  parameters  of  various  controllers  and  the  cost. 

For  ship  speed  15  Knots,  encounter  angle  150.0  degree, 
encounter  frequency  1.50  rad. /sec.  and  sea  state  a. 


CONTR 

K 

Z1 

PI 

P2 

COST 

A 

0.645 

42.40 

07.70 

0 . 0008188 

B  **********  IT  DID  NOx  CONVERGE  ************ 


TABLE  XIV 

Optimal  Controller  Parameters  IX 

Simulation  Results  -  Steady  State  600  seconds 
optimal  parameters  of  various  controllers  and  the  cost. 

For  ship  speed  15  Knots,  encounter  angle  090.0  degree, 
encounter  frequency  0.50  rad. /sec.  and  sea  state  a. 


CONTR  K  Z1  PI  P2  COST 

A  0.53  27.45  08.04  -  0.070697939 


B  **********  IT  Dxd  NOT  CONVERGE  ************ 


TABLE  XV 

Optimal  Controller  Parameters  X 

Simulation  Results  -  Steady  State  600  seconds 
optimal  parameters  of  various  controllers  and  the  cost. 

For  ship  speed  15  Knots,  encounter  angle  090.0  degree, 
encounter  frequency  0.50  rad. /sec.  and  sea  state  6. 


CONTR 

K 

Z1 

El 

P2 

COST 

A 

0.56 

39.71 

11.87 

0.019029424 

B  **********  XT  DID  NOT  CONVERGE  ************ 


TABLE  XVI 

Optimal  Controller  Parameters  XI 


Simulation  Results  -  Steady  State  600  seconds 
optimal  parameters  of  various  controllers  and  the  cost. 

For  ship  speed  10  Knots,  encounter  angle  060.0  degree, 
encounter  frequency  2.50  rad. /sec.  and  sea  state  9. 


CONTR 

K 

Z1 

PI 

P2 

COST 

A 

0.80 

45.00 

10.66 

0.000045194 

B 

0.89 

36.27 

03.22 

03.22 

0.000036077 

TABLE  XVII 

Optimal  Controller  Parameters  XII 

Simulation  Results  -  Steady  State  600  seconds 
optimal  parameters  of  various  controllers  and  the  cost. 

For  ship  speed  20  Knots,  encounter  angle  060.0  degree, 
encounter  frequency  2.50  rad. /sec.  ana  sea  state  9. 


CONTR 

K 

Z1 

PI 

P2 

COST 

A 

0.47 

34.67 

10.35 

0.000014376 

B 

0.61 

23.16 

03.23 

03.03 

0.000011085 

As  seen  from  the  above  results,  for  all  cases  the 
compensators  have  characteristics  of  a  lead  network. 

The  effects  of  sea  state  on  the  controller  parameters 
can  be  seen  by  comparing  the  Tables  XV,  XI  and  XIV  As  seen 
from  the  tables  an  increase  in  sea  state  causes  an  increase 
in  the  cost  value  and  a  decrease  in  the  controller  parame¬ 
ters  as  expected,  because  heavy  sea  state  brings  high 
disturbance  forces  and  moments  and  they  cause  heavy  yawing 
motions.  To  answer  this,  the  controller  time  constants 
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The  results  showed  that  the  cost  curve  is  flat  around 
the  minimum.  The  same  conclusion  is  made  in  [Ref.  7  and  8] 
for  the  SL/7.  Using  controller  parameters  on  the  flat 
portion  of  the  cost  curve,  but  not  at  the  optimum  did  not 
make  a  significant  differences  in  simulation  results,  after 
600  seconds  of  cruising  there  was  little  change  in  the  final 
location  of  the  ship.  This  property  of  the  cost  surface  may 
be  useful  in  reducing  the  time  required  to  find  a  pratical 
minimum  for  the  cost  function. 

The  ship  speed  effects  on  compensator  parameters  for 
regular  seas  can  be  seen  by  comparing  Tables  X,  XVI  and  XVII. 
Observe  that  for  both  cases  increasing  the  ship  speed  made 
the  parameter  values  and  the  cost  value  decrease.  The  same 
observation  can  be  made  by  comparing  the  results  that  are 
found  for  the  SL/7  container  ship  for  different  speeds  in 
[Ref.  9].  These  results  strongly  indicate  that  the  dynamics 
of  the  ship  determines  the  optimum  structure  for  the 
controller . 

To  see  the  stability  of  the  system,  the  ship's  third 
order  Nomoto  model  was  cascaded  with  the  compensator  and  the 
open  loop  system  BODE  plot  was  drawn,  in  every  case  the 
system  is  stable.  The  phase  margin  varies  between  40 
degrees  and  70  degrees  and  the  zero  cross  over  of  the  magni¬ 
tude  curve  is  around  0.04  radian/second.  It  was  observed 
that  small  changes  in  the  compensator  parameters  do  not 
affect  stability.  With  in  large  limits  of  parameter  varia¬ 
tion  the  system  is  always  stable.  For  regular  seas,  ship 
speed  15  knots,  encounter  angle  30.0  degrees,  encounter 
frequency  0.6  radian/ second  and  sea  state  6,  using  compen¬ 
sator  'A'  and  'B' ,  the  structure  of  the  system  is  presented 
in  Figure  (6.12  )  and  the  system  open  loop  BODE  plot  and 
NICHOLS  plot  are  shown  in  Figures  (  6.13  )  and  (  6.14  ). 
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Figure  6.11  The  Cost  Curves  vs.  Z1  and  Pi  when  Parameters 
Changing  Around  the  Optimal  Values 
for  the  Sea  State  Case. 


in 


The  K  values  of  the  curves,  from  top  to  bottom  is 
K  =  0.2 
K  =  0.5 
K  =  0.9 

The  optimal  parameter  values  to  cost  value 
K  =  0 . 53  ,  Zl  =  5 1 . 46  ,  Pl=18.08  and  J=0. 0186974 


Figure  6.10  The  Cost  Curves  vs.  Zl  and  PI  when  Parameters  are 
Changing  Around  the  Optimal  Values 
for  the  Calm  Water  Case. 
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TABLE  XX 

Test  for  LAMBDA  Value  II 


Simulation  Results  -  Steady  State  600  seconds 

for  different  weighting  factor  values, 
regular  seas,  ship  speed  15  knots 
encounter  angle  90.0  degree,  sea  state  8, 
and  encounter  frequency  0.50  radian/ second 
optimal  parameters  of  A  type  controller. 


LAMBDA 

K 

Z1 

PI 

5.41 

0 

71 

24.91 

08.41 

4.91 

0 

68 

25 . 36 

08.37 

4.41 

0 

64 

25.88 

08.32 

3.91 

0 

60 

26.34 

08.22 

3.41 

0 

56 

26.85 

08.11 

2.91 

0 

53 

27.45 

08.04 

2.41 

0 

48 

28.28 

08.04 

1.91 

0 

44 

29.16 

08.16 

1.41 

0 

40 

31.21 

08.35 

0.91 

0 

34 

33.95 

08.71 

0.41 

0 

.26 

39.00 

09.27 

Using  above  compensator  values  did  not  make  a  signifi¬ 
cant  change  on  the  ship  location  after  600  seconds  simula¬ 
tion,  so  it  can  be  said  that  the  accuracy  of  the  weighting 
factor  is  not  very  important. 

A  few  simulation  runs  were  performed  by  changing  the 
optimal  compensator  parameters  a  small  amount  and  the  cost 
curve  was  plotted.  Figure  (  6.10)  shows  the  cost  curves  for 
three  different  K  values  versus  Z1  and  PI  when  the  ship  is 
in  calm  water,  ship  speed  15  knots  and  1  degree  course 
change,  and  Figure  (  6.11  )  shows  the  cost  curves  for  ship 
speed  15  knots,  encounter  angle  60.0  degree,  encounter 
frequency  2.5  radian/ second ,  and  sea  state  9. 


•-  -  — 
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To  see  the  effects  of  the  weighting  factor  (  LAMBDA  ), 
on  the  compensator  parameters,  different  LAMBDA  values,  were 
used  with  compensator  'A'  and  the  compensator  parameters 
were  computed  for  ship  speed  15  knots,  calm  water  and 
regular  seas  with  encounter  angle  90.0  degree,  encounter 
frequency  0.5  radian/ second  and  sea  state  8.  Results  are 
shown  in  Table  (  XIX)  and  Table  (  XX  ) . 


TABLE  XIX 

Test  for  LAMBDA  Value  I 

Simulation  Results  -  Steady  State  600  seconds 

for  different  weighting  factor  values, 
calm  water,  ship  speed  15  knots, 
optimal  parameters  of  A  type  controller. 


LAMBDA 

K 

Z1 

PI 

10.00 

1. 

07 

29.59 

11.30 

9.00 

1. 

00 

31.5 

11.90 

8.00 

0 

.94 

32.91 

12.45 

7.00 

0 

.88 

34.84 

13.14 

6.00 

0 

.80 

37.61 

14.02 

5.41 

0 

.  76 

39.38 

14.62 

4.91 

0 

.72 

41.26 

15 . 19 

4.41 

0 

.68 

43 . 14 

15.78 

3.91 

0 

.63 

45.26 

16.36 

3.41 

0 

.58 

48.06 

17.10 

2.91 

0 

.53 

51.07 

17.93 

2.41 

0 

.48 

55.21 

18 . 9r 

1.91 

0 

.42 

59.79 

20.15 

1.41 

0 

.  35 

64.79 

21.31 

0.91 

0 

.27 

76.89 

24.17 

0.41 

0 

.  18 

89.61 

28.22 
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YOG  (FT) 


Figure  6.9  Simulation  Results  -  Steady  State  600  sec, 


Equation  (  6.1  )  was  added  into  the  fortran  program  and 
the  ship's  motion  were  observed  for  600  seconds  in  calm 
water  and  in  regular  seas  for  30.0  degree  encounter  angle, 
sea  state  6,  encounter  frequency  0.60  radian/ second  and  ship 
speed  15  knots.  Both  compensator  'A'  and  'B'  results  are 
shown  in  Figure  (  6.9  ),  and  after  600  seconds  ship’s  coor¬ 
dinates  are  given  in  Table  XVIII. 

TABLE  XVIII 
Location  of  the  Ship 

Simulation  Results  -  Steady  State  600  seconds 

CASE  Xog  (ft)  Yog  (ft) 

Calm  water 

No  rudder  15200.994873  0.00 

Regular  seas 

Comp. 'A'  15199.0237114  3.2596 

Regular  seas 

Comp.’B’  15199.0763219  2.87116 


For  calm  water  and  regular  seas  compensator  'B'  gave 
better  results  than  compensator  'A'  did,  but  for  some  cases 
compensator  'B'  did  not  converge  and  the  difference  between 
results  is  not  significant.  The  comparisons  were  made  as  to 
which  type  of  compensator  brings  the  ship  nearest  to  final 
location  at  the  end  of  the  600  seconds. 


Figure  6.8  Orientation  of  Space  Axes  and  Moving  Axes 


The  transformation  from  ship  to  space  coordinates  is 
defined  as 

Xog  =  U*cos  (YAW )  -  V -fsin  (YAW)  (eqn  6.1) 

Yog  =  U*sin( YAW)  -  V*cos(YAW) 

Where  Xog,  Yog  =  Coordinates  of  the  center  of  mass  of 

the  ship  relative  to  coordinate  system 
fix  d  with  respect  to  the  surface  of  the 
earth 


46 


As  seen  from  the  figures  the  optimal  parameter  values 
have  a  linear  behavior  when  changing  the  sea  state.  This 
property  would  be  useful  to  create  look  up  tables  for  param¬ 
eter  values.  The  use  of  look  up  tables  will  be  discussed  in 
Chapter  7. 

To  see  the  motion  of  the  ship,  the  ship's  equations  are 
solved  in  ship  coordinates  (x  ,  y)  and  then  transformed  to 
space  coordinates  (x  ,  y),  Figure  (  6.8  )  shows  the  orienta¬ 
tion  of  space  axes  and  moving  axes. 
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To  have  a  better  understanding  about  the  effects  of  sea 
state  on  optimal  controller  parameters  the  cost  value  versus 
parameter  values  were  plotted  for  ship  speed  15  knots, 
encounter  angle  90.0  degree,  encounter  frequency  0.50 
radian/ second  and  sea  state  6  between  8  and  they  are  shown 
in  Figure  (  6.4  ),  (  6.5),  (  6.6)  and  (  6.7). 


YAWE  (DEGREE) 

5  0.50  0.75  1  1-25  1.50 


TIME  (SEC.) 


Figure  6.3  Heading  Errors  in  Degrees  for  Sea  State  6  and  8. 

It  can  be  seen  from  Figure  (  6.2  )  and  (  6.3  ),  for  sea 
state  8  rudder  and  heading  error  values  are  bigger  than  the 
sea  state  6  rudder  and  heading  error  values,  but  it  was 
noticed  that  the  periods  are  almost  the  same  for  both  sea 
states . 
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decrease  and  new  parameters  introduce  more  rudder  motion,  so 
an  increase  in  rudder  motions  and  heading  error  causes 
increase  in  the  cost.  For  ship  speed  15  knots,  encounter 
angle  90.0  degree,  encounter  frequency  0.5  radian/ second , 
sea  state  6  and  sea  state  8,  rudder  angles  and  heading 
errors  were  plotted  in  200  seconds  and  they  are  shown  in 
Figure  (  6.2  )  and  (  6.3). 
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COMP  A 


Where : 


Comp . ' A' 

K - 0*f  34 

Zl---  75.56 
PI---  42.88 
P2 - 


Comp . ’ B 1 
0^30 
61.97 
9.85 
9.72 


Nomoto 
PK=0. 189 
PZ1=  18 . 34 
PP1=7 . 67 
PP2= 116 . 93 


OPEN  LOOP  FREQUENCY  RESPONSE-M08 
USING  ‘8'  TYPE  CASCADE  COMPENSATOR. 


OPEN  LOOP  FREQUENCY  RESPONSE-PHASE 
3  USING  'B'  TYPE  CASCADE  COMPENSATOR. 


Changing  the  environmental  conditions  changes  the 
optimal  controller  parameters.  A  few  simulations  were  run 
changing  the  environmental  conditions  only  but  the  optimal 
controller  parameters  were  kept  unchanged  The  behavior  of 
rudder  and  heading  error  were  observed.  Three  cases  were 
defined  depending  on  the  enviromental  conditions. 

Case  I  ‘.Encounter  angle  90  degree,  encounter  frequency 
0.5  radian/ second ,  sea  state  6  and  ship  speed  15  knots. 
Optimal  controller  parameters  for  Case  I  are 


K  =  0.56 


Z1  =  39.71  PI  =  11.87 


Case  II  : Encounter  angle  90  degree,  encounter  frequency 
0.5  radian/ second ,  sea  state  8  and  ship  speed  15  knots. 
Optimal  controller  parameters  for  Case  II  are 

K  =  0.53  Zl  =  27.45  Pi  =  08.04 

Case  III : Encounter  angle  30  degree , encounter  frequency 
0.5  radian/ second ,  sea  state  6  and  ship  speed  15  knots. 
Optimal  controller  parameters  for  Case  III  are 

K  =  0.358  Zl  =  66.60  PI  =  24.61 

Figure  (  6 . 15  )  and  figure  (  6.16  )  show  rudder  and 
heading  error  when  the  ship  is  in  Case  I  condition  using 
Case  I  and  Case  II  parameters. 

Figure  (  6 . 17  )  and  figure  (  6.18  )  show  rudder  and 
heading  error  when  the  ship  is  in  Case  II  condition  using 
Case  II  and  Case  I  parameters. 

Figure  (  6 . 19  )  and  figure  (  6.20  )  show  rudder  and 
heading  error  when  the  ship  is  in  Case  III  condition  using 
Case  III  and  Case  I  parameters. 

Figure  (  6.21  )  and  figure  (  6.22  )  show  rudder  and 
heading  error  when  the  ship  is  in  Case  I  condition  using 
Case  I  and  Case  III  parameters. 


D  (DEGREE) 

0.50  -0.25  0.00  0.25  0.50  0.75 


Figure  6.15  Rudder  Motion  in  Case  I  with  Case  I 
and  Case  II  Parameters. 
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Figure  6.16  Heading  Error  in  Case  I  with  Case  I 
and  Case  II  Parameters. 


D  (DEGREE) 


o 


Figure  6.18  Heading  Error  in  Case  II  with  Case  II 
and  Case  I  Parameters. 


Figure  6.19  Rudder  Motion  in  Case  III  with  Case  III 

and  Case  I  Parameters. 
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'igure  6.20  Heading  Error  in  Case  III  with  Case  III 
and  Case  I  Parameters. 


Figure  6.21  Rudder  Motion  in  Case  I  with  Case  I 
and  Case  III  Parameters. 


61 


Figure  6.22  Heading  Error  in  Case  I  with  Case  I 
and  Case  III  Parameters. 


As  seen  from  Figures  (  6 . 15  ) ,  (  6.16  ),  (  6 . 17  )  and 
(  6.18  )  using  Case  I  optimal  parameters  in  Case  II  or  Case 
II  optimal  parameters  in  Case  I  did  not  make  a  big  differ¬ 
ence  in  rudder  motion  and  heading  error,  except  in  the  tran¬ 
sient  response  part. 


Figures  (  6.19  )  and  (6.20  )  show  that  to  use  Case  I 


optimal  parameters  when  ship  is  in  Case  III  is  not  proper 
because  those  parameters  increase  rudder  and  heading  error. 


As  seen  from  Figure  (  6.21  ),  using  Case  III  parameters 
in  Case  I  provided  better  rudder  motion  after  the  transient 
response  part.  Cost  values  were  calculated  and  it  was 
observed  that  when  the  ship  is  in  Case  I  using  Case  I 
optimum  parameters  it  gave  a  better  cost  value  than  Case  III 
parameters  did,  end  of  the  600  seconds  simulation.  But  it 
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was  also  observed  that  if  the  final  time  of  simulation  is 
extended,  the  difference  between  the  cost  values  decreases 
and  if  we  continue  to  increase  the  final  time,  using  Case 
III  parameters  in  Case  I  gives  a  smaller  cost  value  than 
Case  I  parameters  do. 

These  results  show  that  the  transient  response  part  is 
important  in  finding  the  optimum  control  parameters. 


VII.  AN  APPROACH  TO  AN  ADAPTIVE  AUTOPILOT 


As  seen  from  the  previous  chapter  the  optimal  controller 
parameters  change  for  changes  in  ship  conditions  such  as 
ship  speed,  loading  etc.  and  also  changes  in  environmental 
conditions  such  as  sea  state,  encounter  angle  encounter 
frequency  and  depth  of  water.  To  maintain  optimal  steering 
performance  automatically  in  the  presence  of  changing  condi¬ 
tions,  it  is  necessary  to  design  an  adaptive  system  which  is 
capable  of  self  adjustment  of  the  controller  parameters  to 
provide  minimum  added  drag  due  to  steering. 

The  steering  control  system,  if  desired  as  an  adaptive 
autopilot,  would  consist  of  four  Subsystems  as  shown  in 
Figure  (  7.1  ) . 

•  Subsystem  #1  would  be  a  computer  which  will  perform  to 
find  the  optimal  control  parameters.  It  should  get  the 
information  about  ship  steering  characteristics  from 
system  state  sensors  such  as  a  gyrocompass,  rudder 
angle  potentiometer  and  speed  log  and  it  should  feed 
the  control  parameter  values  to  the  controller. 

•  Subsystem  #2  would  be  the  controller,  which  should  be 
adjustable.  It  gets  the  parameter  values  from 
subsystem  #1  and  sends  the  rudder  command  to  subsystem 
#3  which  steers  the  ship. 

•  Subsystem  #3  is  the  plant  which  includes  the  system 
state  sensors  and  ship  steering  devices  such  as  servo 
motors,  hydraulic  pumps,  rudder,  steering  gear  etc. 

•  Subsystem  #4  is  a  manual  control  option  for  safety 
rules.  If  manual  steering  is  needed  it  cuts  the 
connection  between  controller  and  steering  devices  and 
gives  the  control  to  the  helmsman. In  case  of  computer 
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failure  it  cuts  the  connection  between  computer  and 
controller  and  sends  to  the  controller  parameter  values 
which  are  chosen  by  the  watch  officer.  It  contains 
look  up  tables  which  specify  control  parameters  as 
functions  of  steering  and  environmental  conditions. 


Figure  7.1  Adaptive  Control  Scheme. 


Success  of  the  this  adaptive  autopilot  system  depends  on 
the  computer  program  as  well  as  the  accuracy  of  the  system 
sensors.  The  computer  program  which  was  used  in  this 
research  may  be  used  on  board,  but  the  present  program  mini¬ 
mization  subroutine  needs  a  lot  of  computation  time  and  it 
also  needs  starting  guesses  for  parameters  which  are 
different  for  every  condition.  If  computation  time  is 
reduced  to  a  reasonable  time  and  starting  values  are  made 
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available  in  proper  limits,  the  modified  function  minimiza¬ 
tion  program  would  be  considered  as  an  adequate  program  for 
on  board  purposes.  The  work  on  reducing  the  computation 
time  is  presented  in  [Ref.  7]. 

In  the  future,  ships  could  have  better  measurement  of 
navigation  than  can  be  provided  by  conventional  equipment  on 
board.  For  example,  the  U.S  Navy  is  involved  in  a  program 
to  build  a  system  which  is  called  NAVSTAR/ GLOBAL  POSITION 
SYSTEM  (GPS).  The  system  will  provide  extremely  accurate 
three-dimensional  position  and  velocity  information  to  users 
anywhere  in  the  world.  And  also  another  system  is  called 
NAVY  REMOTE  OCEAN  SENSING  SYSTEM  (NROSS)  will  be  able  to 
determine  wind  velocities  over  th  world's  oceans  with  an 
accuracy  sufficient  to  determine  ocean  surface  waves.  Using 
such  valid  information  the  watch  officer  can  use  look  up 
tables  and  insert  them  into  the  computer,  so  system  opera¬ 
tion  will  be  very  close  to  the  minimum  cost  value  and  the 
function  minimization  program  will  accomplish  the  fine 
tuning  rapidly.  Detailed  information  about  GPS  and  NROSS 
can  be  found  in  [Ref.  10,  11,  12,  13], 


VIII.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  STUDY 


A.  CONCLUSIONS 

The  conclusions  resulting  from  this  research  of  the 
Mariner  class  ship  fuel  consumption  as  it  relate  to  steering 
might  be  listed  as  follows. 

•  A  steering  control  system  would  minimize  propulsion 
losses  due  to  steering  and  maintain  desired  heading 
with  reasonable  heading  error  in  every  ship  and  enviro- 
mental  condition. 

•  The  study  shows  that  the  best  model  to  describe  the 
dynamics  of  the  ship  is  the  Taylor’s  series  expansion, 
which  allows  both  linear  and  nonlinear  terms  in  the 
ship's  equation  of  motions.  Also  the  third  order  ship 
Nomoto  model  is  reasonable  to  use  instead  of  the  ship's 
equation  of  motions.  It  involves  both  the  sway  and  yaw 
equations . 

•  It  is  believed  that  the  cost  function  which  is 
presented  in  Chapter  6  and  is  commonly  used  by  many 
researchers  is  an  adequate  function  for  on  board  use. 
It  has  variables  such  as  heading  error  and  rudder  angle 
which  can  be  easily  measured  on  board,  and  a  weighting 
factor  (LAMBDA)  which  is  also  easy  to  calculate  but 
depends  on  ship  conditions  such  as  ship  speed.  Results 
in  chapter  7  show  that  accuracy  of  the  LAMBDA  value 
does  not  make  significant  changes  in  the  controller. 

•  In  this  study  two  different  types  of  controller  were 
tried  which  have  been  called  controller  'A'  and  'B'. 
The  structure  of  these  controllers  is  shown  in  Chapter 
7.  Controller  'A'  was  determined  to  be  a  best  struc¬ 
ture,  because  in  some  cases  the  adaptive  calculations 
for  controller  'B'  did  not  converge. 


YAW2  =  YAW2  +  Q2-'DELT 
C  COST  FUNCTION 

TDIFF=TDIFF+  (YAW- YAW2 ) **2 


PROGRAM  TO  CALCULATE  OPTIMAL  GAINS  FOR  CONTROLLER 


//TANSAN  JOB  ( 1789 , 035 6) ,’ RESEARCH CLASS=J 
//"MAIN  0RG=NPGVM1. 1789P 
//  EXEC  FRTXCLGP , IMSL=DP ,REGION= 1024K 
/ / FORT . SYSIN  DD  * 

C  IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 
C  OBTAINED  CHANGE  XS(*)  TO  X(*)  AND  DELETE  XU(*),AND 
XL ( * ) . 

DIMENSION  XS (3), XU (3), XL (3) 

XS ( 1 )  =  0 . 6 
XS (2  1  =41 . 58 
XS ( 3 ) = 11 . 33 

C  XS(I)  IS  THE  STARTING  GUESS 

C  XL (I)  IS  THE  LOWER  LIMIT  FOR  THE  I ' TH  VARIABLE 
C  XU  (I)  IS  THE  UPPER  LIMIT  FOR  THE  I ’ TH  VARIABLE 
XL ( 1 ) = . 3 
XU(1)=0. 9 
XL ( 2 1  =  30 . 

XU  ( 2 1  =  5  0  . 

XL ( 3 )  =  6  . 

XU|  3  1  =  15  . 

C  A  DESCRIPTION  OF  THE  FOLLOWING  PARAMETERS 
C  IS  DISCUSSED  IN  BOXPLX 
R=  9 . / 13  . 

NTA= 1000 
NPR=0 
NAV  =  0 
NV=  3 
IP  =  0 

C  THE  FOLLOWING  STATEMENT  MUST  BE  CHANGED  TO 
C  CALL  PLANT (XX ) 

C  IF  ONLY  SIMULATION  IS  WANTED 

CALL  BOXPLX f NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER) 
WRITE  (6,25) 

25  FORMAT (IX,’  OPTIMAL  GAINS ’,/ ) 

DO  30  1=1,3 

■(6 ,46)1 .XS(I) 
iT(iX,  'x{  ,12,  '  )=  ’  ,F14. 7) 


WRITE(6 ,40)1 .XS(I) 

FORMAT (ix, 'x( ' ,12, ' )= ’ ,F14. 7) 

STOP 

END 

SUBROUTINE  PLANT (XX) 

SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
COMMON  TDIFF 
REAL--8  L,L2  ,L3  ,L4,L5  ,L6 

REAL- 8  X , XDOT , Y , YDOT , U . UDOT , V , VDOT . YAW , R, RDOT 
REAL-8  TIME , ETIME , XI ,X^,X3.^4,X5,X6,X7 ,X^ 

REAL-8  Y0,Y1,Y2,Y3 ,y4,Y^ ,Y^ ,Yt ,Y& 

REAL-8  N0,N1,N2,N3,N4,N5,N6,N7,N8 

REAL-8  C1,C2,C3,C4, C5,F1,F2,F3 

REAL'- 8  R0,DELT,S,DU,U1,K.Z1,Z2,P1,P2 

REAL- 8  DYAWE , YAWE , YAWC , I SR, ISE , TDIFF , LAMDA 

REAL-8  S1.S2 ,DS1,6S2 ,D 

REAL-8  MASS , IZ , XG , YVDOT , NVDOT , YR , YRDOT 

REAL- 8  NR , NRDOT , FX , FY , MZ , RXR , RYR , RXI , RYI , MZR , MZI 

REAL- 8  RX , RY . RZ , TX , TY , TZ , WA , WE 

DIMENSION  XX ( 3 } 

INITIAL  CONDITIONS  FOR  INTEGRATION 
SIMULATION  END  TIME  IN  SECONDS 
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TABLE  XXII 

Sea  State  vs  Wave  Height 


Sea  state 

5 

6 

7 

8 
9 


Wave  height 


5.0 
10.0 
17 . 5 
20.0 
25.0 


(ft) 


The  program  is  set  up  to  calculate  the  optimal  gains  for 
controller  A.  It  can  be  modified  to  obtain  optimal  gains  for 
the  rest  of  the  controllers.  After  obtaining  the  optimal 
gains  the  program  must  be  modified  to  do  a  simulation.  The 
program  has  sufficient  comments  for  appropriate  changes. 

This  program  can  be  modified  to  obtain  the  Nomoto  models. 
It  is  referenced  in  Chapter  4.  The  following  need  to  be 
changed . 

C  GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 
K=XX ( 1 ) 

Z1=XX(2) 

P 1  =  XX ( 3 ) 

P2=XX(4) 

C  ERROR  SIGNAL  TO  DRIVE  RUDDER  (YAW  ACTUAL  -  YAW  COMMAND) 
c  FOR  EQUATIONS  OF  MOTION. 

D= YAW  -  YAWC 

C  ERROR  SIGNAL  TO  DRIVE  RUDDER  (YAW  COMMAND  -  YAW  ACTUAL) 

C  FOR  NOMOTO  3RD  ORDER  MODEL. 

D2= YAWC- YAW2 
DQ1= ( D2  -  Q1 ) / PI 
DQ2=((K*((Z*DQ1)+Q1)/P2 
C  INTEGRATION 

Q1 = Q 1 + DQ1*DELT 
Q2  =  Q2  +  DQ2“DELT 


-/a*./- 
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•  The  real  and  imaginary  value  of  force  or  moment  is  read 
from  the  sea  state  program  output  depending  on  ship 
speed,  encounter  frequency  and  encounter  angle. 

•  These  values  are  converted  to  magnitude  and  phase 
values . 

•  Using  the  formula  below,  forces  and  moment  are  created 
and  added  into  the  ship's  equations  of  motion. 

Force  or  Moment  =  WA*MAGNITUDE*cos (WE*TIME  +  PHASE) 

Where : 

WE  =  Encounter  frequency  ( radian/ second ) 

WA  =  Significant  wave  height  (ft) 


The  correspondence  between  sea  state  and  significant  wave 
height  is  indicated  in  Table  XXI  [Ref.  29].  During  this 
research  the  values  that  are  presented  in  Table  XXII  were 
used  as  significant  wave  height  (WA). 


TABLE  XXI 

Sea  State  vs  Range  for  Wave  Height 


( 


State 
ft  scale ) 


Range  for  .wave . height 
(Feet ) 


0 

1 

2 

3 

4 

5 

6 

7 

8 
9 


0. 

1. 

3. 

6. 

9. 

13 


0.0 
0.32 
65-0 . 98 
96-3.28 
28-4.92 
56-8.20 
84-13  .  1 
1-18.2 


18.2-24.6 

23.0-32.9 
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APPENDIX  C 

PROGRAM  TO  CALCULATE  OPTIMAL  GAINS 


Algorithm  used  here  is  showed  in  Figure  (  C.l  ) 


- Controller 
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APPENDIX  B 

DATA  FOR  SEA  STATE  PROGRAM 


The  sea  state  program  which  is  explained  in  [Ref.  4], 
needs  information  about  the  ship  to  calculate  the  ship's 
added  mass,  added  inertia  and  seaway  disturbance  forces  and 
moments.  Information  about  the  Mariner  was  gathered  from 
[Ref.  27,  28]  and  presented  here  in  the  form  which  the  sea 
state  program  needs.  The  current  line  is  drawn  on  next  page 
to  show  the  location  of  the  values  in  the  format 


■*  ■  .  ■  "'■»  my  ■■■  "  J  ■  4  ’  ■.«  1  u  i 


TDIFF= ISE+ ISR 
GO  TO  200 
400  CONTINUE 
STOP 
END 

ENTRY 

END 
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NO  =  (  0 . 00059) * ( R0*L3*S*S ) 
N0=0 . 0 


Nl= 
N2  = 
N3  = 
N4  = 
N4  = 
N5  = 
N6  = 
N7  = 
N8  = 
COMMON 
Cl  = 


0.0555  )*(R0*L3*S ) 
0.345)*  R0*L3/S) 

r\  r\  o  /  /.  t  v-  /  t*  r\  o  -  -  c 


(RO*L3*S) 
‘ ‘ _  *S ) 
NR - MA  S  S  *X6  *  S ) 

'  0*L4/ 


0 . 00264,  __ 

-0 . 0349 )*{R0*L4 


-1. 158)*(R0*L4/S) 
-0.0293 ) ~ (RO*L3*S*S ) 
0 . 00482 )*£R0*L3*S*S ) 
0.1032)~(RO*L3) 
COEFFICIENTS 
0. 177j*(RO*L3) 


C2= (0.327 
C2=  (MASS-Yv: 


®tL3) 


C3=(-0.00I8 )*[kO*L4 ) 
C3= (MASS*XG- YRDOT ) 


0 . 0 1 7  5 ) * ( R0*L5 ) 
IZ-NRDOT) 

0. 00478 1*(R0*L4) 


C4 
C4 

C5=  .  _  .  _T  ..  w 

C5=  (MAS  S  *XG  -  ftVf)OT  ) 

REGULAR  WAVES 
FX= WA*RX*DC0 S  (WE-'TIME  +  TX 
FY=WA*RY*DCOS (WE*TIME+TY 
MZ=WA*RZ*DCOS (WE*TIME+TZ 

IF (TIME . EQ .0.0)  FX=0,0 
IFfDABS(FY) .LT. 0.00000001)  FY=0.0 
IF(DABSfMZ) .LT. 0.00000001)  MZ=0.0 
EQUATIONS  OF  MOTION 

FI  =  X1*DU  +  X2*DU*DU  +  X3*DU*DU*DU  +  X4*V*V 
1  -  X5*R"R  +  X6*D*D  ♦  X7*V*R  +  X8*V*D  +  FX 

Y0  +  Y1*V  +  Y2*V*V*V  +  Y  3  *V*D*D  +  Y4*R 
Y  5  *R*V*V  +  Y6 "D  +  Y7*D*D*D  +  Y8*D*V*V  +  FY 
NO  +  N1*V  +  N2"V*V*V  +  N3*V*D*D  +  N4-R 

XT  C  -'.-tJ  7  VrT  7  J.  XT£  ^  XT  “7  'Jr r> Tx  r\  XT  Q  -'.-tx-A-t  tV.-t  ;  , 


F2  = 


F3  = 


+  N5*R*V*V  ♦  N 6 *D  +  N7*D*D*D  +  N8*D*V*V  +  MZ 


UDOT 

VDOT 

RDOT 


=  FI/Cl 

-  (  rh* P' 


=  f 


m 


C2*C4-C5*C3 


C2*C4-C5*C3 


GO  TO  50 
GO  TO  50 


'C4*F2-C3*F3 
r  9  *  F  ^  -  C  S  V  9 
C  WHEN  TO  PRINTOUT 

IF  (TIME. EQ. 0.0 
IF  (ICOUNT.EQ.4 
GO  TO  300 
C  CONVERT  RADIANS  TO  DEGREES 
50  YAWDEG=  YAW-57.296 
RDEG=R*57 .296 
RDDEG=RDOT*57 . 296 
DDEG=D*57 .296 
YAWC=YAWC"57 .296 
WRITE  (6,100)  TIME ,  U  ,  V  ,  R 
100  FORMAT C  ,1X,F9.2,1A,£9.6,1X,F9.6,IX,F9.6) 
ICOUNT= 1 

C  TEST  IF  WANT  TO  STOP 

300  IF  (TIME . GT . ETIME )  GO  TO  400 
C  INTEGRATION  STEP  SIZE  DELT 
DELT= 1 . 

C  INTEGRATION 

U=U+UDOT*DELT 
V=V+VDOT*DELT 
R=R+RDOT*DELT 
YAW = YAW +  R*DELT 
S1=S1+DS 1*DELT 

... ... ...  I  ?.  j.  §  ?  *.  5  §  *  R  5  I ... 

TIME = TIME +DELT 
I COUNT  = I COUNT  + 1 
ISE= ISE  ♦  LAMDA*YAWE**2 
ISR= ISR  +  D**2 
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nn  nn  n  nn 


FX=0 . 

FY=0 . 

MZ  =  0  . 

RXR= - . 37 126D5 

RXI= . 68406D5 

RYR= - . 39983D6 

RYI= . 2447D6 

MZR= . 296D8 

MZI=- . 17546D8 

RX= (RXR**2+RXI**2 )** . 5 

RY= ( RYR**2 + RYI**2 1 ** . 5 

RZ=  C MZR**2 + MZ I **2 1 ** . 5 

TX=DATAN2 (RXI ,RXR) 

TY=DATAN2 (RYI ,RYI 1 
TZ=DATAN2 (MZI ,  MZR) 

SIGNIFICANT  WAVE  HEIGHT: (SEA  STATE  2) 

WA=25 . 

ENCOUNTER  FREQUENCY : (WHEN  ENCOUNTER  ANGLE  IS  00) 

WE=0 . 75 

ADDED  MASS  AND  ADDED  INERTIA  TERMS: 

MASS= . 14685D+07 
IZ= .23567D+11 
XG= -23.9 
YR= - . 95066D8 
YRDOT= . 12211D8 
NR= - . 13152D11 
NRDOT=- . 72177D10 
YVDOT= - . 72459D6 
NVDOT  =  - .53846D8 

HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  HERE  AS  PARAMETERS 
RO= 1.987  6* . 5 
YAWE=0 . 0 
DS 1=0 . 0 
DS2=0 . 0 
S1=0.0 
S2=0 . 0 

200  CONTINUE 
INPUT  YAW  COMMAND 
YAWC=0 . 0 

IF  (TIME. GE. 0.0)  YAWC= ( 1 . 0 / 57 . 296 ) 

ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL  -  YAW  ORDERED) 
YAWE = YAW  -  YAWC 
S  =  (  (U*U  )  +  (  V * V ) )  * * .  5 
DU=U-U1 

DS1= (YAWE  -  S1)/P1 
D=(SI  +  DSl*Zlj*K 

AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 
X1=(-0.0253)-;TRO*L2*S) 

X2=(0. 00948 1*(R0*L2) 

X3  =  ( -  0 . 0021/  )  ••  (RO*L2/  S ) 

X4= (-0.189) * ( RO*L2 ) 


fAWE=YAW  -  YAWC 

;=(  (U*u)  +  (V“V)  )••"  .5 

)U=U-U1 

)S1= (YAWE  -  S1)/P1 
)=(SI  +  DS1*ZI)*K 
kL  FORCE  HYDRODYNAMIC  C( 
1 1 = ( -  0 . 0  2  5  3 ) * (RO7 L2*S ) 
12= (0.00948 j*(R0*L2) 

[3= (-0 . 0021/ )~(R0~L2/S) 
C4= (-0.189 )* (RO*L2 ) 

C5= (0 . 00379 )" (R0*l4) 

C6=  C -0 . 02 )*(R0*L2*S"S ) 


X8  = ( 0 . 0 19  6 ) " (RO*L2*S ) 

LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 
Y0=(- 0.0008 )*(RO*L2*S*S ) 

Y0=0 . 0 


Y8= (0 . 25 )*( 
MOMENT  ABOUT  Z 


*fRO*L2*S*. 
51*(RO*L2*; 
R0*L2 ) 


R0*L2 ) 

-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 
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APPENDIX  A 
THESIS  FORTRAN 


$  JOB 


REAL*8  L,L2,L3,L4,L5,L6 

REAL*8  X , XDOT , Y , YDOT , U . UDOT , V , VDOT . YAW , R, RDOT 

REAL "8  TIME , ETIME , XI , X2  ,  X3 .44 ,X5 .X6 ,X7 ,xfi 

REAL* 8  Y0,Yl,Y2,Y3,Y4,Y5,Y6,Yt,Y8 

REAL* 8  N0,N1,N2,N3,N4,N5,N6,N7, N8 

REAL* 8  C1,C2.C3.C4.C5,F1,F2 ,F3 

REAL*8  RO.DELT.S .DO.Ul.K. Z1.Z2 .P1.P2 

REAL* 8  DYAWE , YAWE , YAWC , I SR , ISE , TDIFF , LAMDA 

REAL*8  S1.S2,DS1,DS2,D 

REAL*8  MASS , IZ , XG , YVDOT , NVDOT 

REAL *8  YR , YRDOT , NR , NRDOT , FX , FY , MZ 

REAL* 8  RX , RY , RZ , TX , TY , TZ , WA , WE 

REAL* 8  RXR.RYR.RXI ,RYI .MZR.MZI 

INITIAL  CONDITIONS  FOR  INTEGRATION 

SIMULATION  END  TIME  IN  SECONDS 
ETIME=210. 

TIME=0 . 0 
I COUNT =1 

INITIALIZE  THE  COST  FUNCTION 
ISE=0 . 0 
ISR=0 . 0 
TDIFF=0 . 0 
LAMDA=2 . 91 

GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 
K=0 . 5 
Zl=43 . 32 
Pl=9 . 26 

X , XDOT , Y , YDOT  ARE  FIX  COORDINATES  ON  EARTH 
X=0.0 
Y=0 . 0 
XDOT=0 . 0 
YDOT=0 . 0 

U, UDOT. V. VDOT  ARE  FIX  COORDINATES  ON  SHIP 
Fl=6.6 
F2=0 . 0 
F3=0 . 0 
V  =  0.0 
UDOT=0 . 0 
VDOT=0 . 0 
YAW=0 . 0 
YAWE=0 . 0 
R=0 . 0 
RD0T=0 . 0 

ORDERED  SPEED  IN  FEET /SEC 
15. *1.689  FT/SEC= 15  KNOTS 
Ul= 15. *1.689 

AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 
U=U1 

D  =  RUDDER  ANGLE 
D=0. 0/57. 296 
L=520 , 

L2=L*"2 
L3=L*L"L 
L4=L*L 3 
L5=L*L4 
L6=L*L5 

FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
MOMENTS  IN  Z 


FEET/SEC 
■  15  KNOTS 


70 


The  transient  response  has  a  big  affect  on  finding  the 
optimum  parameter  values.  For  future  studies  it  may  be 
better  to  increase  final  time  so  that  the  effect  of  the 
transient  response  would  not  be  significant. 

Recent  studies  on  roll  stabilization  shows  that  using 
rudder  stabilization  is  successful  in  reducing  roll, 
.[Ref.  14,  15].  Adding  the  roll  equation  into  the  ship  model 
and  determining  a  proper  cost  function  for  minimum  roll 
motion  a  controller  could  be  designed  for  roll 
stabilization. 


•  An  adaptive  controller  which  minimizes  propulsion 
losses  due  to  steering  is  needed  when  ship  characteris¬ 
tics  and  environmental  conditions  change. 

•  For  performance  in  fuel  saving,  an  adaptive  controller 
may  be  better  than  the  existing  Universal  Gyropilot 
(UGP )  . 

•  Since,  equations  of  motion  of  surface  ships  differ  only 
in  the  numerical  values  of  the  coefficients,  the  simu¬ 
lation  programs  used  for  the  Mariner  class  ship  would 
be  useable  for  other  type  of  ships,  knowing  their  hull 
characteristics .  The  studies  made  for  the  SL/7 
containership  are  examples.  [Ref.  6,  9,  7,  8]. 

B.  RECOMMENDATIONS  FOR  FUTURE  STUDY 

This  research  does  not  cover  all  ship  and  enviromentel 
conditions,  so  to  get  a  better  understanding  about  optimal 
controller  parameters,  it  is  recommended  that  the  methods  be 
applied  to  find  the  controller  parameters  for  an  expanded 
range  of  operating  conditions. 

This  thesis  investigated  only  course_keeping  with 
emphasis  on  minimizing  rudder  and  yawing  activity  to  reduce 
fuel  consumption.  If  a  track  following  control  or  an  auto¬ 
matic  control  for  replenishment  at  sea  is  desired,  a 
different  cost  function  might  be  needed.  The  nature  of  the 
cost  function  for  such  applications  should  be  studied. 

Irregular  seas  case  were  not  considered  in  this  study. 
For  future  work  it  is  necessary  to  study  the  effect  of 
irregular  seas  on  the  controller  parameters  and  on  the  cost 
function.  It  is  believed  that  compering  the  regular  sea 
results  with  irregular  sea  results  would  give  better  under¬ 
standing  about  ship's  steering  characteristics. 
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c 

c 

c 

c 


c 

c 

c 

c 


c 

c 


c 

c 

c 


COST  FUNCTION 


OT  ARE  FIX  COORDINATES  ON  EARTH 


0 

,0 


IN 


FEET/SEC 
' '  KNOTS 


ETIME=600 . 

TIME  =  0 . 0 
ICOUNT= I 
INITIALIZE  THE 
ISE  =  0 . 0 
ISR=0 .0 
TDIFF  =  0 . 0 
LAMDA=2 .91 

GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 
K=XX (1) 

Z1=XX(2) 

P1=XX( 3 j 
X , XDOT , Y , YD 

x=o!o 

Y  =  0 . 0 
XDOT=0 . 0 
YDOT=0 . 0 

U , UDOT , V , VDOT  ARE  FIX  COORDINATES  ON  SHIP 
Fl=6.6 
F2  =  0 . 0 
F3  =  0 . 0 

V  =  0.0 
UDOT=0 
VD0T=0 
YAW=0 . 0 
YAWE=0 . 0 
R=0 . 0 
RDOT=0 . 0 

ORDERED  SPEED 
15. *1.689  FT/ SEC= 15 
Ul=15. *1.689 

AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 
U=U1 

D  =  RUDDER  ANGLE 
D=0. 0/57 .296 
L=  520  . 

L2=L**2 
L3=L*L*L 
L4=L*L3 
L5=L*L4 
L6=L*L5 
FORCES  IN  X , Y 
MOMENTS  IN  2 
FX=0  . 

FY=0 . 

MZ  =  0  . 

RXR- . 33374D3 
RXI= . 16064D3 
RYR= - . 17254D5 
RYI= . 1435D4 
MZR= . 30976D7 
MZI= . 846  7  2D6 
RX= ( RXR**2  +  RXI **2 
RY= (RYR**2  +  RYI**2 ' 

RZ= (MZR**2+MZI**2 \ 

TX=DATAN2 (RXI ,RXR 
TY=DATAN2 (RYI , RYI 
TZ=DATAN2 (MZI ,MZR 
SIGNIFICANT  WAVE  ' 

WA=25 . 

ENCOUNTER  FREQUENCY: (WHEN  ENCOUNTER  ANGLE  IS  00) 

WE  =  2 . 5 

ADDED  MASS  AND  ADDED  INERTIA  TERMS: 

MASS= . 14685D+07 
IZ= .23567D+11 
XG= -  2  3  .  9 
YR=- . 24965D8 
YRDOT= . 75199D7 
NR=- .47281D10 


DIRECTION  COMPUTED  IN  FORCES 


5 

5 

5 


/EIGHT:  (SEA  STATE  2) 
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NRDOT  =  - .29233D11 
YVDOT= - . 19723D6 
NVDOT= - . 23 146D8 


C  HYDRODYNAMIC  COEFFICIENTS  ARE 
RO=l. 9876*. 5 
YAWE  =  0 . 0 
DS1=0 . 0 
DS2=0 . 0 
S1=0.0 
S2=0 . 0 

200  CONTINUE 
C  INPUT  YAW  COMMAND 
YAWC=0 . 0 


INSERTED  HERE  AS  PARAMETERS 


IF  (TIME. GE. 0.0)  YAWC= ( 1 . 0 / 57 . 296 ) 


ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL  -  YAW  ORDERED) 


YAWE=YAW  -  YAWC 
S= (7u*u)+ (V*V) ) ** . 5 
DU=U-U1 


DS1= (YAWE 
D=(SI  +  DS1 


PI 

K 


C  AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 

r\  r\  o  c  o  \  -sV  /  rk-.'-T  o  &  c  \  v  ' 


XI=  i 
X2  = 
X3=i 
X4  = 
X5  = 
X6=i 
X7  = 
X8  = 


-0.189 
0.0037, 
-0.02)* 


,S1 
*Zlj* 
DRODY 
)*(RO 
. _WRO  — 
21/ ) "(RO*L 
9 ) * [RO“L2 ) 
79  )*(R0"L4 


r 0,0253) *(RO*L2*S) 

l/S) 


0. 00948 1*(RO*L2) 
-0 . 0021/  * 


X 


s) 


_ ,  ’RO*L2*S’ 

'0. 168)*(RO*L3) 

’0.0196)"  (R0*L2*S  ) 

LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS 


C  Y0  = 

-0.0008 

Y0  =  ( 

).o 

Yl= 

'-0.244) 

Y2  = 

-1.702) 

Y3  = 

-0.0008 

C  Y4= 

-0.105) 

Y4=i 

YR-MASS 

Y5  = 

3 . 2  3  )  *  ( 

Y6  = 

0.0586) 

Y7  = 

-0.0097 

Y8  = 

0.251*1 

C  MOMENT 

ABOUT  Z 

( SWAY ) 


*S) 


6  )  *  ( RO "L2 
9  75  )  "'(RO“ 
" (RO"L2 ) 


;S*S ) 


L2*S*S ) 


-0 . 0555 }*(RO*L3*S ) 
0 . 345 )*(RO*L3/ S) 
0.00264)*(RO*L3"S) 
-0.0349 )*(RO*L4*S) 


NR-MASS-'X 


NO  =(0.00059) * (RO*L3*S*S ) 
N0=0 . 0 
Nl=  '  ‘ 

N2  = 

N3  = 

N4  = 

N4= 

N5  = 

N6=  i 
N7 : 

N8=  i 
COMMON 
Cl=  < 


(YAW) 


1 . 158 ) *TRO*L4 / S ) 

-0 . 0293 ) ~ (RO*Lj*S*S) 
0 . 0  0  4 8  2 ) * ( RO*L  3  *  S  * S ) 
0. 1032)*(R0*L3) 


C2  = 


COEFFICIENTS 
(O. 177)*(RO*L3) 
0.327  WrO*L3) 


) 


C2= (MASS  - YVDOT) 


C3= (-0.0018 )*(RO*L4) 

cl  _  /  ua  r  r  '.Vvc  wn  nn'r  \ 


C3= (MASS*XG- YRD0T ) 


C4= 
C4  = 


0 . 0175 )*(RO*L5 ) 
.IZ-NRDOT) 

C5=(-0. 00478 j*fRO*L4) 

C5= (MASS*XG-NVDOT) 
REGULAR  WAVES 
FX  =  WA'-RX  *  DCO S  ( WE*T IME  +  TX ' 
F  Y = WA*RY  "  DCO  S  ( WE  '"TIME  +  T  Y 
MZ=WA*RZ*DCOS (WE*TIME+TZ 


IF (TIME . EQ . 0 . 0 )  FX=0.0 
IF (DABS (FY) .LT. 0.00000001)  FY=0.0 
IF (DABS (MZ) .LT. 0.00000001)  MZ=0.0 
EQUATIONS  OF  MOTION 

FI  =  XI "DU  +  X2*DU*DU  +  X3 "DU "DU "DU  +  X4*V*V 
+  X5*R*R  +  X6*D*D  +  X7*V*R  +  X8*V*D  +  FX 
F2  =  YO  +  Y1*V  +  Y2*V*V*V  +  Y3*V*D*D  +  Y4*R 
+  Y5*R*V*V  +  Y6--D  +  Y7*D*D*D  +  Y8*D*V*V  + 
F3  =  NO  +  N1"V  +  N2"V"V"V  +  N3*V*D*D  +  N4*R 
+  N5*R*V*V  +  N6-D  +  N7*D*D*D  ♦  N8*D*V*V  + 


)/ (C2*C4-C5"C3) 
)/ (C2*C4-C5*C3) 


UDOT  =  Fl/Cl 
VDOT  =  (C4*F2-C3*F3 
RDOT  =  (C2*F3-C5*F2 

WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.21)  GO  TO  50 
GO  TO  300 

CONVERT  RADIANS  TO  DEGREES 
50  YAWDEG=  YAW* 5 7. 29 6 
RDEG=R*57 .296 
RDDEG=RDOT*57 .296 
DDEG=D*57 .296 
YAWC=YAWC*57 . 296 
WRITE  (6,100)  TIME , TDIFF 
100  FORMAT ( '  ' ,1X,F10.2,1X,F20. 10) 
ICOUNT= 1 

TEST  IF  WANT  TO  STOP 
300  IF  ( TIME . GT . ETIME )  GO  TO  400 

INTEGRATION  STEP  SIZE  DELT 
DELT= 1 . 

INTEGRATION 

U=U+UDOT*DELT 

V=V+VDOT*DELT 

R=  R+ RDOT*DELT 

YAW=YAW+R*DELT 

S1=S1+DS1*DELT 

i?=??+DS?*DELf 


FY 

MZ 


TIME= TIME + DELT 
ICOUNT=ICOUNT+ 1 
ISE=ISE  +  L AMD A* Y AWE  *  *  2 
ISR=ISR  +  D**2 
GO  TO  200 
400  TDIFF=ISE+ISR 

WRITE (6 , 500 )  TIME. TDIFF 
500  FORMAT (’  T , IX , F9 . 2 , IX , F20 . 10 ) 

RETURN 

END 

C  DELETE  ALL  THE  FOLLOWING  SUBROUTINE  IF  SIMULATION  ONLY 
C  AND  NOT  OPTIMIZATION  IS  WANTED 

C  . 

c 

C  SUBROUTINE  BOXPLX  (CATEGORY  HO) 

C 

C  PURPOSE 

C 

C  BOXPLX  IS  A  SUBROUTINE  USED  TO  SOLVE  THE  PROBLEM  OF 
C  LOCATING  A  MINIMUM  ?OR  MAXIMUM)  OF  AN  ARBITRARY  OBJECT- 
C  IVE  FUNCTION  SUBJECT  TO  ARBITRARY  EXPLICIT  AND/OR 
C  IMPLICIT  CONSTRAINTS  BY  THE  COMPLEX  METHOD  OF  M.J.  BOX. 

C  EXPLICIT  CONSTRAINTS  ARE  DEFINED  AS  UPPER  AND  LOWER 
C  BOUNDS  ON  THE  INDEPENDENT  VARIABLES  IMPLICIT  CONSTRAINTS 
C  MAY  BE  ARBITRARY  FUNCTION  OF  THE  VARIABLES.  TWO  FUN- 
C  CTION  SUBPROGRAM  TO  EVALUATE  THE  OBJECTIVE  FUNCTION  AND 
C  IMPLICIT  CONSTRAINTS,  RESPECTIVELY,  MUST  BE  SUPPLIED 
C  BY  THE  USER  7 SEE  EXAMPLE  BELOWJ .  BOXPLX  ALSO  HAS  tHE 
C  OPTION  TO  PERFORM  INTEGER  PROGRAMMING,  WHERE  THE  VALUES 
c  OF  THE  INDEPENDENT  VARIABLES  ARE  RESTRICTED  TO  INTEGERS. 
C 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


USAGE 

CALL  BOXPLX  (NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER ) 
DESCRIPTION  OF  PARAMETERS 

NV  AN  INTEGER  INPUT  DEFINING  THE  NUMBER  OF  INDEPENDENT 
VARIABLES  OF  THE  OBJECTIVE  FUNCTION  TO  BE  MINIMIZED. 
NOTE:  MAXIMUM  NV  +  NAV  IS  PRESENTLY  50.  MAXIMIM  NV  IS 

25.  IF  THESE  LIMITS  MUST  BE  EXCEEDED.  PUNCH  A  SOURCE 
DECK  IN  THE  USUAL  MANNER,  AND  CHANGE  THE  DIMENSION 
STATEMENTS . 

NAV  AN  INTEGER  INPUT  DEFINING  THE  NUMBER  OF  AUXILIARY  var 
iaBLES  THE  USER  WISHES  TO  DEFINE  FOR  HIS  OWN  CONVENIENCE. 
TYPICALLY  HE  MAY  WISH  TO  DEFINE  THE  VALUE  OF  EACH  IMPLICI 
CONSTRAINT  FUNCTION  AS  AN  AUXILIARY  VARIABLE.  IF  THIS 
IS  DONE,  THE  OPTIONAL  OUTPUT  FEATURE  OF  BOXPLX  CAN  BE 
USED  TO  OBSERVE  THE  VALUES  OF  THOSE  CONSTRAINTS  AS  THE 
SOLUTION  PROGRESSES.  AUXILIARY  VARIABLES,  IF  USED, 

SHOULD  BE  EVALUATED  IN  FUNCTION  KE  (DEFINED  BELOW). 

NAV  MAY  BE  ZERO. 

NPR  INPUT  INTEGER  CONTROLLING  THE  FREQUENCY  OF  OUTPUT 
desired  for  diagnostic  purposes. 

IF  NPR  .LE.  0,  NO  OUTPUT  WILL  BE 

PRODUCED  BY  BOXPLX.  OTHERWISE,  THE  CURRENT  COMPLEX  OF 
K=  2*NV  VERTICES  AND  THEIR  CENTROID  WILL  BE  OUTPUT  AFTER 
EACH  NPR  PERMISSIBLE  TRIALS.  THE  NUMBER  OF  TOTAL  TRIALS, 
NUMBER  OF  FEASIBLE  TRIALS,  NUMBER  OF  FUNCTION  EVALUATIONS 
AND  NUMBER  OF  IMPLICIT  CONSTRAINT  EVALUATIONS  ARE  IN¬ 
CLUDED  IN  THE  OUTPUT. 

ADDITIONALLY,  (WHEN  NPR  .GT.  0)  THE  SAME  INFORMATION 
WILL  BE  OUTPUT: 

IF  THE  INITIAL  POINT  IS  NOT  FEASIBLE, 

AFTER  THE  FIRST  COMPLETE  COMPLEX  IS  dENERATED^ 

IF  A  FEASIBLE  VERTEX  CANNOT  BE  FOUND  AT  SOME  'i'RIAL , 

IF  THE  OBJECTIVE  VALUE  OF  A  VERTEX  CANNOT  BE  MADE 
NO-LONGER- WORST. 

5)  IF  THE  LIMIT  ON  TRIALS  (NTA)  IS  REACHED  AND, 

6)  WHEN  THE  OBJECTIVE  FUNCTION  HAS  BEEN  UNCHANGED  FOR 
2"NV  TRIALS,  INDICATING  A  LOCAL  MINIMUM  HAS  BEEN 
FOUND . 


IF  THE  USER  WISHES  TO 
A  CHOICE  OF  NPR  =  25, 


TRACE  THE  PROGRESS  OF  A  SOLUTION, 
50  OR  100  IS  RECOMMENDED. 


NTA  INTEGER  INPUT  OF  LIMIT  ON  THE  NUMBER  OF  TRIALS 
allowed  in  the  calculation. 

IF  THE  USER  INPUTS  NTA  . LE .  0,  A  default 
VALUE  OF  2000  IS  USED.  WHEN  THIS  LIMIT  IS  REACHED 
CONTROL  RETURNS  TO  THE  CALLING  PROGRAM  WITH  THE  BEST 
ATTAINED  OBJECTIVE  FUNCTION  VALUE  IN  YMN,  AND  THE  BEST 
ATTAINED  SOLUTION  POINT  IN  XS . 

R  A  REAL  NUMBER  INPUT  TO  DEFINE  THE  FIRST  RANDOM  NUMBER 
USED  IN  DEVELOPING  THE  INITIAL  COMPLEX  OF  2*NV  VERTICIES. 
(0.  .GT.  R  .LT.  1.)  IF  R  IS  NOT  WITHIN  THESE  BOUNDS, 

IT  WILL  BE  REPLACED  BY  1 . / 3 .  . 

XS  INPUT  REAL  ARRAY  DIMENSIONED  AT  LEAST  NV+NAV. 
the  first  .iv  must  contain  a 
FEASIBLE  ORIGIN  FOR  STARTING  THE  CAL¬ 
CULATION.  THE  LAST  NAV  NEED  NOT  BE  INITIALIZED.  UPON 
RETURN  FROM  BOXPLX,  THE  FIRST  NV  ELEMENTS  OF  THE  ARRAY 
CONTAIN  THE  COORDINATES  OF  THE  MINIMUM  OBJECTIVE 
function,  AND  THE  REMAINING  NAV  (NAV  .GE.  0)  CONTAIN  THe 
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C  values  of  THE  CORRESPONDING  AUXILIARY  VARIABLES. 

C 

C  IP  INTEGER  INPUT  FOR  OPTIONAL  INTEGER  PROGRAMMING. 

C  if  ip= 1 ,  THE  VALUES  OF  THE  INDEPENDENT  VARIABLES  WILL 
C  be  replaced  WITH  INTEGER  VALUES  (STILL  STORED  AS  REAL*4). 

C  XU  A  REAL  ARRAY  DIMENSIONED  AT  LEAST  NV  INPUTTING  THE 
C  upper  BOUND  ON  EACH  INDEPENDENT  VARIABLE,  (EACH  EXPLICIT 
C  conSTRAINT).  INPUT  VALUES  ARE  SLIGHTLY  ALTERED  BY  BOXPLX. 


C  conSTRAINT).  INPUT  VALUES  ARE  SLIGHTLY  ALTERED  BY  BOXPLX. 
C 

C  XL  A  REAL  ARRAY  DIMENSIONED  AT  LEAST  NV  INPUTTING  THE 
c  lower  bound  on  each  independent 
C  VARIABLE,  (EACH  EXPLICIT  CONstraint). 

C  NOTE:  FOR  BOTH  XU  AND  XL  CHOOSE  REASONABLE 
C  VALUES  IF  NONE  ARE  GIVEN,  NOT  VALUES  WHICH  ARE 
C  magnitudes  ABOVE  OR  BELOW  THE  EXPECTED  SOLUTION. 

C  input  values  are  SLIGHTLY  ALTERED  BY  BOXPLX. 

C 

C  YMN  THIS  OUTPUT  IS  THE  VALUE  (REAL*4)  OF  THE  OBJECTIVE 
C  funcTION, CORRESPONDING  TO  THE  SOLUTION  POINT  OUTPUT  IN  XS 
C 

C  IER  INTEGER  ERROR  RETURN.  TO  BE  INTERROGATED  UPON 
C  return  FROM  BOXPLX.  IER  WILL  BE  ONE  OF  THE  FOLLOWING: 

C 

C  =-l  CANNOT  FIND  FEASIBLE  VERTEX  OR  FEASIBLE  CENTROID 
C  AT  THE  START  OR  A  RESTART  (SEE  'METHOD'  BELOW). 

C  =0  FUNCTION  VALUE  UNCHANGED  FOR  fN'  TRIALS.  (WHERE 
C  N=  6*NV+ 10 )  THIS  IS  THE  NORMAL  RETURN  PARAMETER. 

C  =1  CANNOT  DEVELOP  FEASIBLE  VERTEX. 

C  =  2  CANNOT  DEVELOP  A  NO-LONGER-WORST  VERTEX. 

C  =3  LIMIT  ON  TRIALS  REACHED.  (NTA  EXCEEDED) 

C  NOTE:  VALID  RESULTS  MAY  BE  RETURNED  IN  ANY  OF  THE 
C  ABOVE  CASES . 


C  =1 
C  =2 
C  =3 
C  NOTE 
C 


JF  THE 


EXAMPLE  OF  USAGE 


C  THIS  EXAMPLE  MINIMIZES  THE  OBJECTIVE  FUNCTION  SHOWN  IN 
C  the  EXTERNAL  FUNCTION  FE (X) .  THERE  ARE  TWO  INDEPENDEN1 
C  varlABLES  X(l)  &  X?2),  AND  TWO  IMPLICIT  CONSTRAINT 
C  function  X(3)  &  X(4)  WHICH  ARE  EVALUATED  AS  AUXILIARY 
C  variables  (see  EXTERNAL  FUNCTION  KE(X)  ). 

C  DIMENSION  XS(4) ,XU(2) ,XL(2) 


STARTING  GUESS 
XS?1)  =  1.0 
XS]2j  =0.5 
UPPER  LIMITS 
XU(1)  =  6.0 

xm2j  =  6.0 

LOWER  LIMITS 
XL ( 1 )  =  0.0 
XL(2)  =  0.0 


R  =  9./13. 
NTA  =  5000 
NPR  =  50 
NAV  =  2 
NV  =  2 
IP  =  0 

CALL  BOXPLX 
WRITE76 . 1) 
1FORMAT  (/// 


3?4^eAND  iHE  Motion 


(NV .NAV , NPR . NTA , R , XS , IP , XU , XL , YMN , IER ) 

(xstn,i=1.4)  ,YtiN  ieA) 

, ,  THE  P6lNf  IS  LOCATED  AT  (XS(I)  =  )  ’ 
VALUE  IS  ' , E13 . 7 , '  IER  =  ’,15) 


STOP 

END 
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cc 

cc 

C  FUNCTION  KE(X) 

C  EVALUATE  CONSTRAINTS.  SET  KE=0  IF  NO  IMPLICIT  CONSTRAINT 
C  is  viOLATED,  OR  SET  KE=1  IF  ANY  IMPLICIT 
c  constraint  is  violated. 

C 
C 
C 
C 
C 
C 
C 
C 

CC 
C 
C 
C 

CC 
CC 
C 
C 

CC 
CC 

c 
c 
c 
u 
c 
c 

c  THE  COMPLEX  METHOD  IS  AN  EXTENSION  AND  ADAPTION  OF 
c  the  simple  method  of  linear  programming. 

C  STARTING  WITH  ANY  ONE  feasible  point  in  n-dimension 
C  A  COMPLEX”  OF  2“N  vertices  is  constructed  by 
C  SELECTING  RANDOM  POINTS  WITHIN  THE  feasible 
C  REGION.  FOR  THIS  PURPOSE  N  COORDINATES  ARE  FIRST 
C  RANDOMLY  CHOSEN  WITHIN  THE  SPACE  BOUNDED  BY  EXPLICIT  CON- 
C  STRAINTS.  THIS  DEFINES  A  TRIAL  INITIAL  VERTEX, 
c  it  is  then  checked  for  possible  violation 
C  OF  IMPLICIT  CONSTRAINTS!  IF  one  or  more  are  violated, 

C  THE  TRIAL  INITIAL  VERTEX  IS  DISPLACED  half  of  its 
C  DISTANCE  FROM  THE  CENTROID  OF  PREVIOUSLY  SELECTED  initial 
C  VERTICES.  IF  NECESSARY  THIS  DISPLACEMENT  PROCESS  IS 
C  REPEATED  UNTIL  THE  VERTEX  HAS  BECOME  FEASIBLE.  IF  THIS 
c  fails  to  happen  after  5"n+10  displacements, 

C  THE  SOLUTION  IS  ABANDoned.  after  each  vertex  is  added 
C  TO  THE  COMPLEX,  THE  CURRENT  centroid  is  checked  for 
C  FEASIBILITY.  IF  IT  IS  INFEASIBLE,  the  last  trail 
C  VERTEX  IS  ABANDONED  AND  AN  EFFORT  TO  GENERATE  an  alter- 
C  ATIVE  TRIAL  VERTEX  IS  MADE.  IF  5*N+10  VERTICES  ARE 
C  ABANDONED  CONSECUTIVELY,  THE  SOLUTION  IS  TERMINATED. 

C 

C  IF  AN  INITIAL  COMPLEX  IS  ESTABLISHED,  THE  BASIC 
c  computation  loop  is  initiated. 

C  THESE  INSTRUCTIONS  FIND  THE  CURRENT  WORST  vertex,  that 
C  IS,  THE  VERTEX  WITH  THE  LARGEST  CORRESPONDING  value  for 
C  THE  OBJECTIVE  FUNCTION,  AND  REPLACE  THAT  VERTEX  BY 
C  ITS  OVER-REFLECTION  THROUGH  THE  CENTROID  OF  ALL  OTHER 
c  vertices,  (if  the  vertex  to  be 
C  REPLACED  IS  CONSIDERED  AS  A  VECTOR  IN  n- space, 

C  ITS  OVER-REFLECTION  IS  OPPOSITE  IN  DIRECTION,  IN- 
C  CREASED  IN  LENGTH  BY  THE  FACTOR  1.3,  AND  COLLINEAR  WITH 
C  THE  REPLACED  VERTEX  AND  CENTROID  OF  ALL  OTHER  VERTICES.) 

c 

C  WHEN  AN  OVER-REFLECTION  IS  NOT  FEASIBLE  OR  REMAINS 
c  worst,  it  is  considered  not-permissible 
C  AND  IS  DISPLACED  HALFWAY  TOWARD  the  centroid. 

C  AFTER  FOUR  SUCH  ATTEMPTS  ARE  MADE  UNSUCCESSFULLY 
C  EVERY  FIFTH  ATTEMPT  IS  MADE  BY  REFLECTING  THE  OFFENDING 


DIMENSION  X(4) 

XI  =  X(l) 

X2  =  X(2) 

KE  =  0 

xm  =  XI  +  1.7  3205 1*X2 
IF  (X(3)  .LT.  0.  .OR.  X ( 3 ) 
X (4 )  =  Xl/1. 732051  -X2 
IF  (X(4)  .GE.  0.)  RETURN 

1  KE  =  1 
RETURN 
END 


.GT.  6. )  GO  TO  1 


FUNCTION  FE (X) 
DIMENSION  X (4 ) 


THIS  IS  THE  OBJECTIVE  FUNCTION. 
FE=-(X(2)**3  *(9.-(X(l)-3.  )** 
RETURN 


2)/(46. 76538)) 


END 
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c  vertex  through  the  present  best 

C  VERTEX,  INSTEAD  OF  THROUGH  THE  CENtroid.  if  5*n+10 
C  DISPLACEMENTS  AND  OVER- REFLECTIONS  OCCUR  without  a 
C  SUCCESSFUL  (PERMISSIBLE)  RESULT,  THE  CURRENT  BEST 
C  VERTEX  IS  TAKEN  AS  AN  INITIAL  FEASIBLE  POINT  FOR  A 
c  restart  run  of  the  complete  process. 

C  RESTARTING  IS  ALSO  UNDERTAKEN  when  6*nv+10  consecutive 
C  TRIALS  HAVE  BEEN  MADE  WITH  NO  SIGNIFicant  change  in  the 
C  VALUE  OF  THE  OBJECTIVE  FUNCTION.  IN  ALL  cases, 

C  RESTARTING  IS  INHIBITED  IF  THE  LAST  RESTART  DID  NOT 
C  pRODUCE  A  SIGNIFICANT  IMPROVEMENT  IN  THE  MINIMUM 
c  attained. 

C 

C  IT  IS  RECOMMENDED  THAT  THE  USER  READ  THE  REFERENCE  FOR 
C  FURTHER  USEFUL  INFORMATION.  IT  SHOULD  BE  NOTED  THAT  THE 
C  ALGORITHM  DEFINED  THERE  HAS  BEEN  ALTERED  TO  FIND  THE 
C  CONSTRAINED  MINIMUM,  RATHER  THAN  THE  MAXIMUM. 

C 

C 

C 

C  REMARKS 
C 

C  THE  INTEGER  PROGRAMMING  OPTION  WAS  ADDED  TO  THIS  PROGRAM 
C  AS  SUGGESTED  IN  REFERENCE  (2).  A  MIXED 
c  integer/ continuous  variable  version  of  boxplx 
C  WOULD  BE  EASY  TO  CREATE  BY  DEclaring  nip"  to  be  an  array 
C  OF  NV  CONTROL  VARIABLES  WHERE  IP  7i7=l  would  indicate 
C  THAT  THE  I-TH  VARIABLE  IS  TO  BE  CONFINED  to  integer 
C  VALUES.  EACH  STATEMENT  OF  THE  FORM  'IF  (IP  .EQ.IV  etc. 

C  WOULD  THEN  NEED  TO  BE  ALTERED  TO  'IF  (IP  Cl)  .EQ.  1)'  etc. 
C  ,  WHERE  THE  SUBSCRIPT  IS  APPROPRIATELY  CHOSEN.  NORMALLY, 
C  XU  AND  XL  VALUES  ARE  ALTERED  TO  BE  AN  EPSILON  'WITHIN' 
c  actual  values 

C  DECLARED  BY  THE  USER.  THIS  ADJUSTMENT  IS  NOT  MADE 
C  WHEN  IP=1. 

C 

C  NOTE:  NO  NON-LINEAR  PROGRAMMING  ALGORITHM  CAN  GUARANTEE 
c  that  the  answer  found  is  the  global 

C  MINIMUM,  RATHER  THAN  JUST  A  local  minimum,  however, 

C  ACCORDING  TO  REF. 2.  THE  COMPLEX  method  has  an  advantage 
C  IN  THAT  IT  TENDS  TO  FIND  THE  GLOBAL  minimum  more 
C  FREQUENTLY  THAN  MANY  OTHER  NON-LINEAR  PROGRAM- 
C  MING  ALGORITHMS. 

C 

C  IT  SHOULD  BE  NOTED  THAT  THE  AUXILIARY  VARIABLE  FEATURE 
c  can  also  be  used  to  deal  with 

C  PROBLEMS  CONTAINING  EQUALITY  CONstraints.  any  equality 
C  CONSTRAINT  IMPLIES  THAT  A  GIVEN  VARiable  is  not  truly 
C  INDEPENDENT.  THEREFORE,  IN  GENERAL,  ONE  variable 
C  INVOLVED  IN  AN  EQUALITY  CONSTRAINT  CAN  BE  RENUMBERED  from 
C  THE  SET  OF  NV  INDEPENDENT  VARIABLES  AND  ADDED  TO  THE  SET 
C  OF  NAV  AUXILIARY  VARIABLES.  THIS  USUALLY  INVOLVES 
C  renumbering  THE  INDEPENDENT  VARIABLES  OF  THE  GIVEN 
C  problem 

C  SUBROUTINES  AND  FUNCTIONS  REQUIRED 

C  SUBROUTINE  'BOUT'  AND  FUNCTION  ' FBV '  ARE  INTEGRAL 
C  parts  of  THE  BOXPLX  PACKAGE. 

C 

C  TWO  FUNCTIONS  MUST  BE  SUPPLIED  BY  THE  USER.  THE  FIRST, 
c  ke(x),  is  used  to  evaluate  the  implicit 

C  CONSTRAINTS.  SET  KE=0  AT  THE  beginning  of  the  function 
C  .  THEN  EVALUATE  THE  IMPLICIT  CONstraints.  in  the  example 
C  ABOVE,  THE  FIRST  CONSTRAINT,  X(3J.  must  be  within  the 
C  RANGE  (0.  .LE.  X(3)  .LE.  6.J.  THE  SECOND  constraint  x(4) 
C  ,  MUST  BE  .GE.  0.  .  IF  EITHER  CONSTRAINT  IS  not  within 
C  THESE  BOUNDS,  CONTROL  IS  TRANSFERRED  TO  STATEMENT  1, 

C  AND  KE  IS  SET  TO  nl"  AND  CONTROL  IS  RETURNED  TO  BOXPLX. 

C 
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THE  SECOND  FUNCTION  THE  USER  MUST  PROVIDE  EVALUATES  THE 
objective  function,  it  is 

CALLED  FE (X)  AS  SHOWN  IN  THE  EXAMple  above,  and  fe 
MUST  BE  SET  TO  THE  VALUE  OF  THE  OBJECTIVE  function 
CORRESPONDING  TO  CURRENT  VALUES  OF  THE  NV  INDEPENDENT 
VARIABLES  IN  ARRAY  ' XT . 
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SUBROUTINE  BOXPLX  (NV ,NAV , NPR , NTZ , RZ , XS , IP , BU , BL , YMN , IER) 


DIMENSION  Vf 50,50 ) ,  FUN(50),  SUM(25),  CEN(25),  XS(NV) 
l,bu(nv) ,bl(NV) 


KV  =  5 
EP  =  l.E-6 
NTA  =  2000 
IF  (NTZ.GT.O)  NTA  = 

IF~(R.LE.O. .OR.R.GE.l. )  R=l./3. 
NVT  =  NV  +  NAV 


NTZ 


TOTAL  VARS,  EXPLICIT  PLUS  IMPLICIT 

=  0 

CURRENT  TRIAL  NO. 

=  0 

CURRENT  NO.  OF  PERMISSIBLE  TRIALS 
NTFS  =  0 

CURRENT  NO.  OF  TIMES  F  HAS  BEEN  ALMOST  UNCHANGED 


NT 

NPT 


CHECK  FEASIBILITY  OF  START  POINT 


DO 

VT 

IF 

II 

VT 


4  1=1, NV 


=  XS 
(BL| 


A.  •  LI 


)  .LE.VT)  GO  TO  1 


=  BL(I) 

GO  TO  2 

IF  (BU(I) .GE.VT)  GO  TO  3 
II  =  I 
VT  =  BU(I) 

IF  (NPR. GT . 0 )  WRITE  (6,49) 

** ' '  VT 
VT 

, . 1 )  GO  TO  4 

SET! ) +AMAX1 (EP , EP*ABS (BL(I ) ) ) 

BUm-AMAXl(EP>EP*ABS(BU(I)  )  ) 


II 


NCE  =  1 

NUMBER  OF  CONSTRAINT  EVALUATIONS 

IF~ (KE(V( 1 , 1) ) . EQ. 0 )  GO  TO  5 
IF  (NPR.LE.O)  GO  TO  12 
WRITE  (6,50) 

GO  TO  12 
5  NFE  =  1 


NUMBER  OF  VERTICES  (K)  =  2  TIMES  NO.  OF  VARIABLES. 

K  =  2“NV 

NUMBER  OF  DISPLACEMENTS  ALLOWED. 

NLIM  =  5*NV+ 10 

NUMBER  OF  CONSECUTIVE  TRIALS  WITH  UNCHANGED  FE  TO 
terminate . 

NOT  =  NLIM+NV 
ALPHA  =  1.3 
FK  =  K 
FKM  =  FK-1. 

BETA  =  ALPHA* 1. 

INSURE  SEED  OF  R..NDOM  NUMBER  GENERATOR  IS  ODD. 

I OR  =  R " 1 . E  7 

IF  (MOD ( IQR , 2 ) . EQ . 0 )  IQR=IQR+L01 

SET  UP  INITIAL  VERTICES 
FUN ( 1 )  =  FE (V ( 1 , 1 ) ) 

YMN  =  FUN ( 1 ) 

6  FI  =  1. 

FUNOLD  =  FUN ( 1 ) 

DO  15  1=2, K 
FI  =  FI+1. 

LIMT  =  0 

7  LIMT  =  LIMT+ 1 

END  CALCULATION  IF  FEASIBLE  CENTROID  CANNOT  BE  FOUND. 
IF  ( LIMT. GE. NLIM)  GO  TO  11 

DO  8  J  = 1 , NV 

RANDOM  NUMBER  GENERATOR  (RANDU) 

IQR  =  IQR*65539 

IF  (IQR.LT.O)  IQR  =  IQR+2147483647+ 1 
RQX  =  IQR 

RQX  =  RQX* .4656613E-9 

VCJ.I)  =  BL( Jj+RQX*(BU( J ) -BL(J ) ) 

IF  (iP.EQ.l)  V(J,I)=AINT(V(J,I)+ .5) 

8  CONTINUE 


DO  10  L= 1 , NLIM 
NCE  =  NCE+1 

IF  (KE (V ( 1 , I ) ) . EQ . 0 )  GO  TO  13 
DO  9  J  =  1 . NV 

VT  =  .5*(V(J,I)+CEN(J) ) 

IF  (IP.EQ.I)  VT  =  AlNT (VT+ . 5 ) 


9  cfellJ 


=  VT 
UE 


10  CONTINUE 


11  IF  (NPR.LE.O)  GO  TO  12 
WRITE  (6,5 1)  I 

CALL  BOUT  (NT,NPT,NFE,NCE,NV,NVT,V,I , 


12  IER  =  -1 
GO  TO  48 


FUN , CEN , I  ) 
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13  DO  14  J  =  1 ,  N V 

SUM(J)  =  SUM ( J ) +V ( J , I ) 

14  CEN(J)  =  SUM(J)/FI 

TRY  TO  ASSURE  FEASIBLE  CENTROID  FOR  STARTING. 
NCE  =  NCE+ 1 

IF  (KE(CEN) .EQ.O)  GO  TO  60 
SUMiJ)  =  SUM(J)  -V ( J , I ) 

GO  TO  7 

60  NFE  =  NFE+ 1 


C 

C 

c 


FUN (I )  =  FE ( V( 1 , I ) ) 
15  CONTINUE 


END  OF  LOOP  SETTING  OF  INITIAL  COMPLEX. 

IF  (NPR.LE.O)  GO  TO  17 

CALL  BOUT  (NT , NPT , NFE , NCE , NV , NVT ,  V  ,K , FUN , CEN , 0 ) 

FIND  THE  WORST  VERTEX,  THE  'J’TH. 

J  =  1 

DO  16  1=2, K 

IF  ( FUN ( J ) . GE . FUN ( I ) )  GO  TO  16 

16  CONTINUE 

BASIC  LOOP.  ELIMINATE  EACH  WORST  VERTEX  IN  TURN, 
it  must  become  NO  LONGER  WORST,  NOT  MERELY  IMPROVED, 
find  next- to-vertex,  THE  'JN'TH  ONE. 

17  JN  =  1 

IF  (J.EQ.l)  JN  =  2 


DO  18  1=1, 
IF  (I.EQ.J 
IF  (FUN ( JN 
JN  =  I 
18  CONTINUE 


K 


)  GO  TO  18 

J . GE . FUN ( I ) )  GO  TO  18 


LIMT  =  NUMBER  OF  MOVES  DURING  THIS  TRIAL  TOWARD  THE 
centroid  DUE  TO  FUNCTION  VALUE. 

LIMT  =  1 

COMPUTE  CENTROID  AND  OVER  REFLECT  WORST  VERTEX. 


DO  19  1=1, NV 
VT  =  V(I,J) 

SUM  ( I )  =  SUM  (I )  -1 


CEN  (I 

VT  =  _  __ 

IF  (IP.EQ.l) 


. . , . .  VT 

=  SUM (I ) / FKM 
ETA"'CEN^IJ  -  ALPHA  "'VT 


AINT (VT  +  .  5  ) 


INSURE  THE  EXPLICIT  CONSTRAINTS  ARE  OBSERVED. 

19  V(I,J)  =  AMAX1 ( AMIN1 (VT , BU ( I ) ) , BL ( I ) ) 

NT  =  NT+ 1 

CHECK  FOR  IMPLICIT  CONSTRAINT  VIOLATION. 

20  DO  25  N= 1 , NLIM 
NCE  =  NCE+ 1 

IF  (KE(V(1, J) ) .EQ.O)  GO  TO  26 

EVERY  'KV’TH  TIME,  OVER- REFLECT  THE  OFFENDING  VERTEX 
through  the  BEST  VERTEX. 

IF  (MOD(N.KVl.NE.O)  GO  TO  22 
CALL  FBV  (K,FUN,M) 

DO  21  1=1, NV 

VT  =  BETA" V  ( I ,  M )  -  ALPHA"'V  ( I ,  J  ) 
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on  o  noon  no  '  on  ono  n  n  non 


IF  (IP.EQ.I)  VT  =  AINT ( VT  +  .  5  } 

21  V(I,J)  =  AMAX1 (AMIN1 (VT , BU ( I  )  )  , 

GO  TO  24 

CONSTRAINT  VIOLATION: 


BL(I)  ) 


MOVE  NEW  POINT  TOWARD  CENTROID. 


22  DO  23  1=1. NV 

VT  =  . 5*(CEN(I)+V(I , J 

I?-U?-EQ^1)  vf  -A" 


V 

23  C 


VT+  .  5  ) 


UE 


24  NT  =  NT+ 1 

25  CONTINUE 

IER  =  1 

CANNOT  GET  FEASIBLE  VERTEX  BY  MOVING  TOWARD  CENTROID, 

OR  BY  OVER-REFLECTING  THRU  THE  BEST  VERTEX. 

IF  (NPR.LE.O)  GO  TO  42 
WRITE  (6,52)  NT , J 

CALL  BOUT  (NT , NPT ,NFE , NCE ,NV , NVT , V ,K , FUN , CEN , J ) 

GO  TO  42 

FEASIBLE  VERTEX  FOUND,  EVALUATE  THE  OBJECTIVE  FUNCTION. 

26  NFE  =  NFE+1 
FUNTRY  =  FE ( V ( 1 , J ) ) 

TEST  TO  SEE  IF  FUNCTION  VALUE  HAS  NOT  CHANGED. 

AFO  =  ABS  (FUNTRY- FUNOLD) 

AMX  =  AMAXl(ABS(EP*FUNOLD) ,EP) 

ACTIVATE  THE  FOLLOWING  TWO  STATEMENTS  FOR  DIAGNOSTIC 
purposes  only. 

WRITE  (6,99) 

AFO .AMX .FUNTRY .FUNOLD , FUN ( J ) , FUN ( JN ) , NTFS , N 
9$  FORMAT  (l$,I3,6Ei5.7j2l$) 

IF  (AFO.GT.AM&)  GO  t6  27 
NTFS  =  NTFS+ 1 
IF  (NTFS.LT.NCT)  GO  TO  28 
IER  =  0 

IF  (NPR.LE.O)  GO  TO  42 
WRITE  (6,53)  K 
CALL  BOUT  (NT,1 
GO  TO  42 

27  NTFS  =  0 


,  NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , 0 ) 


C 

C 

C 

C 


IS  THE  NEW  VERTEX  NO  LONGER  WORST? 
28  IF  (FUNTRY. LT. FUN ( JN ) )  GO  TO  34 


TRIAL  VERTEX  IS  STILL  WORST:  ADJUST  TOWARD  CENTROID. 
EVERY  'KV  TH  TIME,  OVER-REFtECT  THE  OFFENDING 


C 

c 


hi  f  kliVA  IN  T  All  11.11U  t  UYU1N 

through  the  BEST  VERTEX 
LIMT  =  LIMT+ 1 
IF  (MOD(LIMT.KV) .NE.O)  GO  TO  30 
CALL  FBV  (K,FUN,M) 

DO  29  1=1, NV 

VT  =  BETA"V ( I , M) - ALPHA "V (I , J ) 

IF  (IP.EQ.l)  $T  =  AINT ( VT + . 5  j  , 

29  V(I,J)  =  AMAX1 ( AMIN1 ( VT , BU ( I ; ) , BL ( I ) ) 

GO  TO  32 

30  DO  31  1=1, NV 

VT  =  .  5*(CEN(I) J V(I , J) ) 

IF  (IP.EQ.l)  VT  =  AiNT (VT+ . 5 ) 

V(I,J)  =  VT 


VERTEX 


31  CONTINUE 

32  IF  (LIMT.LT.NLIM)  GO  TO  33 

CANNOT  MAKE  THE  ’ J ’ TH  VERTEX  NO  LONGER  WORST  BY 
displacing  toward 

THE  CENTROID  OR  BY  OVER- REFLECTING  THRU  THE  BEST  VERTEX 
IER  =  2 

IF  (NPR  .LE.  0)  GO  TO  42 
WRITE  (6 ,52 j  NT,  J 

CALL  BOUT  (NT,NPf , NFE , NCE , NV , NVT , V ,K , FUN , CEN , J) 

GO  TO  42 

33  NT  =  NT+ 1 
GO.  TO  20 

SUCCESS:  WE  HAVE  A  REPLACEMENT  FOR  VERTEX  J. 

34  FUNTjj  =  FUNTRY 
FUNOLD  =  FUNTRY 
NPT  =  NPT+ 1 

EVERY  100 'TH  PERMISSIBLE  TRIAL ,  RECOMPUTE  CENTROID 
summation  to  AVOID  CREEPING  ERROR. 

IF  (MOD (NPT , 100 ) . NE . 0 )  GO  TO  37 

DO  36  1=1. NV 
SUM ( I )  =  0. 

DO  35  N=1,K 

35  SUM(I)  =  SUM (I)+V(I,N) 

CEN (I)  =  SUM ( I ) / FK 

36  CONTINUE 

LC  =  0 
GO  TO  39 

37  DO  38  1=1, NV 

38  SUM ( I )  =  0UM(I)+V(I , J) 

LC  =  J 


39  IF  (NPR.LE.O)  GO  TO  40 

IF  (MOD (NPT , NPR) . NE . 0 )  GO  TO  40 

CALL  BOUT  ( NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , LC ) 

HAS  THE  MAX.  NUMBER  OF  TRIALS  BEEN  REACHED  WITHOUT 
convergence?  iF  NOT,  GO  TO  NEW  TRIAL. 

40  IF  (NT.GE.NTA)  GO  TO  41 

NEXT-TO-WORST  VERTEX  NOW  BECOMES  WORST. 

J  =  JN 
GO  TO  17 

41  IER  =  3 

IF  (NPR.GT.O)  WRITE  (6,54) 

COLLECTOR  POINT  FOR  ALL  ENDINGS. 

1)  CANNOT  DEVELOP  FEASIBLE  VERTEX.  IER  =  1 

2)  CANNOT  DEVELOP  A  NO -LONGER- WORST  VERTEX.  IER  =  2 

3)  FUNCTION  VALUE  UNCHANGED  FOR  K  TRIALS.  IER  =  0 

4)  LIMIT  ON  TRIALS  REACHED.  IER  =  3 

5j  CANNOT  FIND  FEASIBLE  VERTEX  AT  START.  IER  =  -1 

42  CONTINUE 

FIND  BEST  VERTEX. 

CALL  FBV  (K.FUN.M) 

IF  (IER.GE.3)  GO  TO  44 

RESTART  IF  THIS  SOLUTION  IS  SIGNIFICANTLY  BETTER  THAN 


c 

c 


the  previous,  OR  IF  THIS  IS  THE  FIRST  TRY. 

IF  (NPR.LE.O)  GO  TO  43 
WRITE  (6,55)  (M, YMN , FUN(M) ) 

43  IF  ( FUN (Mj . GE . YMN )  GO  TO  47 

IF  (ABS (FUN(M) - YMN ) . LE . AMAX1 (EP , EP*YMN) )  GO  TO  47 


GIVE  IT  ANOTHER  TRY  UNLESS  LIMIT  ON  TRIALS  REACHED. 
44  YMN  =  FUN?Mj 
FUN ( 1 )  =  FUN (M) 


DO  45  1=1, NV 
CEN(I)  =  v(i,m; 
SUM  £  1 1  =  V(I,M 
45  V ( I , 1 )  =  V  I,M 


DO  46  1=1, NVT 
46  XS(I)  =  V(l,M) 


IF  (IER.LT.3)  GO  TO  6 

47  IF  (NPR.LE.Ol  GO  TO  48 

CALL  BOUT  ( NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , V ( 1 , M ) 
WRITE  (6,56)  £UN(A) 

48  RETURN 


1) 


49  FORMAT  (50H0INDEX  AND  DIRECTION  OF  OUTLYING 
lvariable  at  starti5) 

50  FORMAT  (50H0IMPLICIT  CONSTRAINT  VIOLATED  AT 
lstart.  dead  end.) 

51  FORMAT  ('OCANNOT  FIND  FEASIBLE 14 ,’ TH  VERTEX  OR 
iid  r*"  _J - *-  ' 


lcentroid  at  start.  ) 

52  FORMAT  ( 10H0AT  TRIAL  I4,54H  CANNOT  FIND  FEASIBLE 
lvertex  which  is  no 

2LONGER  W0RST.I4, 15X, ’RESTART  FROM  BEST  VERTEX.’) 

53  FORMAT  (40HO^UNfiTI01^  HAS  BEEN  ALMOST  UNCHANGED 
lfor  i5 , 7h  trails] 

54  FORMAT  (27H0LIMIT  ON  TRIALS  EXCEEDED.  ) 

55  FORMAT (  uBEST  VERTEX  IS  NO.’,I3,’OLD  MIN  WAS’,E15.7, 
1  T  NEW  MIN  IS  ’ , E15 .7) 

56  FORMAT  ( ’ OMIN  OBJECTIVE  FUNCTION 
END 

SUBROUTINE  FBV  (K , FUN ,M) 

DIMENSION  FUN (50) 

M  =  1 


IS  ’ , E15 . 7 ) 


DO  1  1=2, K 

IF_ (FUN(M) .LE.FUN(I) )  GO  TO 
CONTINUE 


RETURN 

END 

SUBROUTINE  BOUT  (NT , NPT , NFE , NCE ,NV ,NVT , V , K , FN , C , IK ) 
DIMENSION  V  ( 50 , 5u]  ,  FN(.^O),  C(2^) 

WRITE  (6,4)  NT, NPT, NFE, NCE 


C 

C 


ffItNV&l?NV)(Gi4o(l’I),J=1’NV) 

NVP  =  NV+ 1 

WRITE  (6,6)  (V(J,I) , J=NVP,NVT) 
CONTINUE 


IF  (IK.NE.O)  GO  TO  2 


WRITE  (6,7)  (C(I) ,1=1, NV ) 
RETURN 

IF  (IK.GE.O)  GO  TO  3 
WRltE  (6,8)  (C(I) , I = 1 , NV ) 
RETURN 


\ 
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3  WRITE  (6,9)  IK, (C(I) ,1=1, NV ) 

RETURN 

C 

4  FORMAT  ('ONO.  TOTAL  TRIALS  =  ',I5,4X, 
l’no.  feasible  trails  =  ’ ,i5,4x, 

2 ' NO .  FUNCTION  EVALUATIONS  =  ',I5,4X, 

3'no.  constraint  evaluations  =  ',i5 / 

4 ' 0  FUNCTION  VALUE '  6X  ' INDEPENDENT  VARIABLES/ 

5dependENT  OR  IMPLICIT  CONSTRAINTS') 

5  FORMAT  (1H  ,  E18 . 7 , 2X , 7E 14 . 7 / ( 21X , 7E14 . 7 ) ) 

6  FORMAT  (21X.7E14J) 

7  FORMAT  ( lOHOCENTROID  11X , 7E14 . 7 / ( 2 IX . 7E14 . 7 ) ) 

8  FORMAT  rO  BEST  VERTEX  ’  7X  ,  7E14 . 7  /  (  2 IX  ,  7E14 . 7  )  ) 

9  FORMAT  f'OCENTROID  LESS  ^X’  ,I2.,2X,  7E14. 7/ (21X, 
17el4 . 7  )  ) 

END 

FUNCTION  FE (X ) 

DIMENSION  X(3) 

COMMON  TDIFF 
CALL  PLANT (X) 

FE=TDIFF 

RETURN 

END 

FUNCTION  KE (X ) 

DIMENSION  X(3) 

KE  =  0 

RETURN 

END 

/ l GO . SYSIN  DD  * 

L  THIS  CARD  SHOULD  BE  USED  ONLY  REGULAR  SEA  CASE 
/ / GO . FT12F00 1  DD  DISP= SHR , DSN=MSS . S2160 . A341 
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